Meta-Analyzing ICC Estimates Graphs

rm(list = ls())
library(tidyverse)
library(kableExtra)
library(ggplot2)


load("data/results_pooled_est.RData")

load("data/results_var_comp.RData")

1 Table

table <- perf_crit  |>
  group_by(model, variance, var_combo_graph) |>
  summarize(
    mean_pb = mean(bias_icc),
    var_pb = var(bias_icc),
    min_pb = min(bias_icc),
    max_pb = max(bias_icc),
    
    # mean_pb_mcse = mean(bias_icc_mcse),
    # var_pb_mcse = var(bias_icc_mcse),
    # min_pb_mcse = min(bias_icc_mcse),
    # max_pb_mcse = max(bias_icc_mcse),
    
    
    mean_rpb = mean(rel_bias_icc),
    var_rpb = sd(rel_bias_icc),
    min_rpb = min(rel_bias_icc),
    max_rpb = max(rel_bias_icc),
    
    mean_rmse = mean(RMSE),
    var_rmse = sd(RMSE),
    min_rmse = min(RMSE),
    max_rmse = max(RMSE),
    
    mean_mse = mean(MSE),
    var_mse = sd(MSE),
    min_mse = min(MSE),
    max_mse = max(MSE),
    
    mean_rseb = mean(rel_bias_se_ICC),
    var_rseb = sd(rel_bias_se_ICC),
    min_rseb = min(rel_bias_se_ICC),
    max_rseb = max(rel_bias_se_ICC),
    
    # mean_rpb_mcse = mean(rel_bias_icc_mcse),
    # var_rpb_mcse = var(rel_bias_icc_mcse),
    # min_rpb_mcse = min(rel_bias_icc_mcse),
    # max_rpb_mcse = max(rel_bias_icc_mcse
    .groups = "drop"
  ) |>
  arrange(var_combo_graph) |>
  rename(
    Estimation = model,
    "Variance Formulae" = variance,
    "ICC Parameter" = var_combo_graph,
    "Parameter Bias - mean" = mean_pb,
    "Parameter Bias - variance" = var_pb,
    "Parameter Bias - min." = min_pb,
    "Parameter Bias - max" = max_pb,
    "Relative Parameter Bias - mean" = mean_rpb,
    "Relative Parameter Bias - variance" = var_rpb,
    "Relative Parameter Bias - min." = min_rpb,
    "Relative Parameter Bias - max" = max_rpb,
    "Root Mean Squared Error - mean" = mean_rmse,
    "Root Mean Squared Error - variance" = var_rmse,
    "Root Mean Squared Error - min." = min_rmse,
    "Root Mean Squared Error - max" = max_rmse,
    "Mean Squared Error - mean" = mean_mse,
    "Mean Squared Error - variance" = var_mse,
    "Mean Squared Error - min." = min_mse,
    "Mean Squared Error - max" = max_mse,
    "Relative Standard Error Bias - mean" = mean_rseb,
    "Relative Standard Error Bias - variance" = var_rseb,
    "Relative Standard Error Bias - min." = min_rseb,
    "Relative Standard Error Bias - max" = max_rseb#,
    
    #
    # mean_rpb_ratio = mean(rel_bias_icc_alt),
    # min_rpb_ratio = min(rel_bias_icc_alt),
    # max_rpb_ratio = max(rel_bias_icc_alt),
    # sd_rpb_ratio = sd(rel_bias_icc_alt)
    
  )

  write.csv(table, file = "perform_crit_summary.csv")
  
  table |> 
   kbl() |>
   kable_styling() |> 
   scroll_box(width = "900px", height = "500px")

Estimation	Variance Formulae	ICC Parameter	Parameter Bias - mean	Parameter Bias - variance	Parameter Bias - min.	Parameter Bias - max	Relative Parameter Bias - mean	Relative Parameter Bias - variance	Relative Parameter Bias - min.	Relative Parameter Bias - max	Root Mean Squared Error - mean	Root Mean Squared Error - variance	Root Mean Squared Error - min.	Root Mean Squared Error - max	Mean Squared Error - mean	Mean Squared Error - variance	Mean Squared Error - min.	Mean Squared Error - max	Relative Standard Error Bias - mean	Relative Standard Error Bias - variance	Relative Standard Error Bias - min.	Relative Standard Error Bias - max
REML	Donner	0.05	-0.0055206	3.10e-06	-0.0086502	-0.0033179	-0.1104120	0.0350042	-0.1730048	-0.0663589	0.0068103	0.0018541	0.0038674	0.0109862	0.0000497	0.0000262	1.50e-05	0.0001207	-0.0418569	0.0298111	-0.0980127	0.0093520
REML	Fisher	0.05	-0.0055134	3.10e-06	-0.0084706	-0.0032662	-0.1102687	0.0351356	-0.1694117	-0.0653237	0.0068032	0.0018560	0.0038551	0.0108288	0.0000497	0.0000261	1.49e-05	0.0001173	-0.0419004	0.0297949	-0.0973853	0.0089611
REML	Fisher TF	0.05	-0.0038014	1.30e-06	-0.0055870	-0.0020606	-0.0760285	0.0224458	-0.1117401	-0.0412123	0.0053695	0.0015421	0.0026996	0.0083844	0.0000312	0.0000170	7.30e-06	0.0000703	-0.0021009	0.0235978	-0.0595003	0.0393137
REML	Hedges	0.05	-0.0054133	3.00e-06	-0.0084060	-0.0032763	-0.1082661	0.0343600	-0.1681204	-0.0655256	0.0067228	0.0018241	0.0038339	0.0108028	0.0000485	0.0000253	1.47e-05	0.0001167	-0.0406881	0.0295053	-0.0949330	0.0108128
REML	Smith	0.05	-0.0056378	3.20e-06	-0.0088097	-0.0033974	-0.1127560	0.0358532	-0.1761949	-0.0679483	0.0069050	0.0018834	0.0039056	0.0110766	0.0000512	0.0000269	1.53e-05	0.0001227	-0.0428401	0.0300361	-0.0980784	0.0085310
REML	Swiger	0.05	-0.0055107	3.10e-06	-0.0084615	-0.0032660	-0.1102139	0.0350982	-0.1692307	-0.0653193	0.0068010	0.0018545	0.0038549	0.0108251	0.0000496	0.0000261	1.49e-05	0.0001172	-0.0418829	0.0297847	-0.0971925	0.0089368
RVE	Donner	0.05	-0.0055777	3.00e-06	-0.0085649	-0.0033146	-0.1115539	0.0346341	-0.1712981	-0.0662929	0.0068763	0.0018542	0.0038117	0.0109456	0.0000507	0.0000262	1.45e-05	0.0001198	-0.0541107	0.0453614	-0.1367120	0.0244123
RVE	Fisher	0.05	-0.0055726	3.00e-06	-0.0083736	-0.0033254	-0.1114524	0.0347894	-0.1674713	-0.0665089	0.0068696	0.0018546	0.0037992	0.0107733	0.0000506	0.0000261	1.44e-05	0.0001161	-0.0532168	0.0450808	-0.1369035	0.0266975
RVE	Fisher TF	0.05	-0.0038014	1.30e-06	-0.0055868	-0.0020606	-0.0760289	0.0224448	-0.1117352	-0.0412113	0.0053696	0.0015422	0.0026996	0.0083850	0.0000312	0.0000170	7.30e-06	0.0000703	-0.0020099	0.0235566	-0.0594670	0.0397569
RVE	Hedges	0.05	-0.0054668	2.90e-06	-0.0083203	-0.0032798	-0.1093365	0.0339513	-0.1664052	-0.0655957	0.0067853	0.0018239	0.0037795	0.0107628	0.0000493	0.0000254	1.43e-05	0.0001158	-0.0540002	0.0443273	-0.1361919	0.0245538
RVE	Smith	0.05	-0.0056991	3.20e-06	-0.0087006	-0.0033549	-0.1139811	0.0355106	-0.1740112	-0.0670980	0.0069730	0.0018801	0.0038499	0.0110027	0.0000521	0.0000269	1.48e-05	0.0001211	-0.0530666	0.0462642	-0.1369892	0.0268108
RVE	Swiger	0.05	-0.0055698	3.00e-06	-0.0083648	-0.0033253	-0.1113955	0.0347507	-0.1672966	-0.0665057	0.0068672	0.0018531	0.0037991	0.0107691	0.0000505	0.0000261	1.44e-05	0.0001160	-0.0532047	0.0450506	-0.1368946	0.0266430
REML	Donner	0.10	-0.0070980	4.50e-06	-0.0103775	-0.0043971	-0.0709803	0.0211630	-0.1037751	-0.0439705	0.0086478	0.0022452	0.0050989	0.0130228	0.0000797	0.0000395	2.60e-05	0.0001696	-0.0277921	0.0230186	-0.0795846	0.0331137
REML	Fisher	0.10	-0.0071406	4.60e-06	-0.0103323	-0.0044362	-0.0714056	0.0214930	-0.1033234	-0.0443619	0.0086834	0.0022679	0.0051412	0.0131500	0.0000804	0.0000401	2.64e-05	0.0001729	-0.0281511	0.0229956	-0.0802938	0.0325047
REML	Fisher TF	0.10	-0.0049637	1.70e-06	-0.0073567	-0.0025374	-0.0496370	0.0130568	-0.0735668	-0.0253741	0.0068912	0.0017562	0.0034347	0.0109561	0.0000505	0.0000250	1.18e-05	0.0001200	-0.0026267	0.0205246	-0.0413956	0.0444442
REML	Hedges	0.10	-0.0069785	4.40e-06	-0.0101173	-0.0043557	-0.0697847	0.0209693	-0.1011728	-0.0435567	0.0085501	0.0022251	0.0050570	0.0128375	0.0000780	0.0000386	2.56e-05	0.0001648	-0.0270140	0.0231271	-0.0794918	0.0341317
REML	Smith	0.10	-0.0073212	4.80e-06	-0.0108172	-0.0044442	-0.0732116	0.0220040	-0.1081715	-0.0444417	0.0088344	0.0023074	0.0051495	0.0131570	0.0000833	0.0000415	2.65e-05	0.0001731	-0.0288332	0.0233585	-0.0802933	0.0325404
REML	Swiger	0.10	-0.0071367	4.60e-06	-0.0103242	-0.0044352	-0.0713670	0.0214699	-0.1032416	-0.0443520	0.0086802	0.0022661	0.0051401	0.0131430	0.0000804	0.0000400	2.64e-05	0.0001727	-0.0281266	0.0229980	-0.0802703	0.0325331
RVE	Donner	0.10	-0.0070986	4.50e-06	-0.0102447	-0.0043907	-0.0709858	0.0212475	-0.1024472	-0.0439069	0.0086500	0.0022356	0.0051176	0.0129352	0.0000797	0.0000393	2.62e-05	0.0001673	-0.0126710	0.0331201	-0.0803805	0.0558477
RVE	Fisher	0.10	-0.0071480	4.70e-06	-0.0102092	-0.0044334	-0.0714798	0.0216354	-0.1020923	-0.0443341	0.0086913	0.0022620	0.0051666	0.0130718	0.0000805	0.0000398	2.67e-05	0.0001709	-0.0123514	0.0335987	-0.0793201	0.0568207
RVE	Fisher TF	0.10	-0.0049635	1.70e-06	-0.0073566	-0.0025376	-0.0496354	0.0130562	-0.0735663	-0.0253758	0.0068911	0.0017563	0.0034347	0.0109561	0.0000505	0.0000250	1.18e-05	0.0001200	-0.0030384	0.0205266	-0.0416797	0.0444772
RVE	Hedges	0.10	-0.0069702	4.40e-06	-0.0099708	-0.0043473	-0.0697021	0.0209913	-0.0997078	-0.0434734	0.0085450	0.0022116	0.0050710	0.0127432	0.0000778	0.0000382	2.57e-05	0.0001624	-0.0137747	0.0322466	-0.0801877	0.0542139
RVE	Smith	0.10	-0.0073410	4.90e-06	-0.0107213	-0.0044420	-0.0734098	0.0222217	-0.1072128	-0.0444202	0.0088525	0.0023053	0.0051759	0.0130782	0.0000836	0.0000413	2.68e-05	0.0001710	-0.0106762	0.0344995	-0.0792582	0.0655423
RVE	Swiger	0.10	-0.0071438	4.70e-06	-0.0102015	-0.0044323	-0.0714385	0.0216096	-0.1020149	-0.0443229	0.0086878	0.0022600	0.0051653	0.0130641	0.0000805	0.0000398	2.67e-05	0.0001707	-0.0123816	0.0335700	-0.0793435	0.0567584
REML	Donner	0.15	-0.0081131	5.70e-06	-0.0114168	-0.0051556	-0.0540873	0.0159411	-0.0761119	-0.0343704	0.0100361	0.0025978	0.0061445	0.0152725	0.0001073	0.0000531	3.78e-05	0.0002332	-0.0244915	0.0239332	-0.0795617	0.0473734
REML	Fisher	0.15	-0.0081880	5.90e-06	-0.0116082	-0.0051958	-0.0545870	0.0162508	-0.0773881	-0.0346387	0.0100986	0.0026305	0.0061904	0.0152935	0.0001088	0.0000540	3.83e-05	0.0002339	-0.0248300	0.0240200	-0.0803382	0.0469676
REML	Fisher TF	0.15	-0.0063898	3.70e-06	-0.0095806	-0.0032772	-0.0425987	0.0127480	-0.0638704	-0.0218479	0.0086216	0.0022393	0.0043784	0.0137339	0.0000792	0.0000399	1.92e-05	0.0001886	-0.0077120	0.0244446	-0.0627136	0.0337777
REML	Hedges	0.15	-0.0079890	5.70e-06	-0.0111627	-0.0051144	-0.0532603	0.0158501	-0.0744181	-0.0340960	0.0099354	0.0025843	0.0060984	0.0150839	0.0001053	0.0000521	3.72e-05	0.0002275	-0.0237531	0.0239181	-0.0797239	0.0472985
REML	Smith	0.15	-0.0083992	6.20e-06	-0.0119814	-0.0052049	-0.0559948	0.0166630	-0.0798761	-0.0346993	0.0102747	0.0026759	0.0062003	0.0156848	0.0001126	0.0000559	3.84e-05	0.0002460	-0.0255192	0.0242548	-0.0841158	0.0469529
REML	Swiger	0.15	-0.0081836	5.90e-06	-0.0115976	-0.0051948	-0.0545571	0.0162335	-0.0773176	-0.0346319	0.0100949	0.0026284	0.0061892	0.0152860	0.0001087	0.0000539	3.83e-05	0.0002337	-0.0248006	0.0240173	-0.0803106	0.0469787
RVE	Donner	0.15	-0.0081877	6.00e-06	-0.0115896	-0.0051975	-0.0545850	0.0163212	-0.0772643	-0.0346503	0.0100978	0.0026211	0.0062413	0.0152677	0.0001087	0.0000537	3.90e-05	0.0002331	-0.0006777	0.0294194	-0.0746190	0.0730476
RVE	Fisher	0.15	-0.0082726	6.30e-06	-0.0116808	-0.0052429	-0.0551507	0.0166864	-0.0778721	-0.0349526	0.0101690	0.0026597	0.0062968	0.0152887	0.0001103	0.0000547	3.97e-05	0.0002337	-0.0001122	0.0298060	-0.0747036	0.0738986
RVE	Fisher TF	0.15	-0.0063894	3.70e-06	-0.0095800	-0.0032776	-0.0425958	0.0127463	-0.0638666	-0.0218505	0.0086213	0.0022391	0.0043784	0.0137329	0.0000792	0.0000399	1.92e-05	0.0001886	-0.0082602	0.0244407	-0.0634913	0.0330649
RVE	Hedges	0.15	-0.0080501	5.90e-06	-0.0112830	-0.0051529	-0.0536675	0.0161842	-0.0752201	-0.0343530	0.0099856	0.0026029	0.0061880	0.0150658	0.0001063	0.0000526	3.83e-05	0.0002270	-0.0017014	0.0290064	-0.0748050	0.0720901
RVE	Smith	0.15	-0.0085070	6.70e-06	-0.0122454	-0.0052528	-0.0567136	0.0171966	-0.0816357	-0.0350185	0.0103658	0.0027158	0.0063082	0.0157074	0.0001147	0.0000570	3.98e-05	0.0002467	0.0017423	0.0305753	-0.0730971	0.0740775
RVE	Swiger	0.15	-0.0082676	6.20e-06	-0.0116683	-0.0052417	-0.0551172	0.0166663	-0.0777885	-0.0349448	0.0101648	0.0026574	0.0062954	0.0152810	0.0001102	0.0000547	3.96e-05	0.0002335	-0.0001479	0.0297885	-0.0747047	0.0738751
REML	Donner	0.25	-0.0086569	6.50e-06	-0.0117173	-0.0055425	-0.0346275	0.0102083	-0.0468692	-0.0221701	0.0114385	0.0029811	0.0068325	0.0168380	0.0001395	0.0000693	4.67e-05	0.0002835	-0.0133702	0.0258710	-0.0560011	0.0644447
REML	Fisher	0.25	-0.0087611	6.80e-06	-0.0118725	-0.0055795	-0.0350445	0.0104341	-0.0474900	-0.0223181	0.0115223	0.0030168	0.0068786	0.0169508	0.0001417	0.0000705	4.73e-05	0.0002873	-0.0137947	0.0258242	-0.0560912	0.0642463
REML	Fisher TF	0.25	-0.0089397	9.20e-06	-0.0134210	-0.0044222	-0.0357587	0.0121512	-0.0536839	-0.0176887	0.0116512	0.0031704	0.0057559	0.0177719	0.0001456	0.0000743	3.31e-05	0.0003158	0.0010399	0.0253765	-0.0712062	0.0528747
REML	Hedges	0.25	-0.0085434	6.50e-06	-0.0114737	-0.0054110	-0.0341736	0.0101869	-0.0458949	-0.0216439	0.0113503	0.0029770	0.0067886	0.0166920	0.0001375	0.0000686	4.61e-05	0.0002786	-0.0127728	0.0258112	-0.0561434	0.0640485
REML	Smith	0.25	-0.0089778	7.10e-06	-0.0124868	-0.0055883	-0.0359112	0.0106921	-0.0499473	-0.0223533	0.0116977	0.0030550	0.0071475	0.0172391	0.0001460	0.0000724	5.11e-05	0.0002972	-0.0146388	0.0259807	-0.0562466	0.0648744
REML	Swiger	0.25	-0.0087564	6.80e-06	-0.0118605	-0.0055775	-0.0350257	0.0104238	-0.0474420	-0.0223099	0.0115185	0.0030150	0.0068773	0.0169447	0.0001416	0.0000705	4.73e-05	0.0002871	-0.0137662	0.0258247	-0.0560896	0.0642467
RVE	Donner	0.25	-0.0088005	7.00e-06	-0.0119670	-0.0055898	-0.0352018	0.0105707	-0.0478682	-0.0223592	0.0115555	0.0030182	0.0069429	0.0168908	0.0001425	0.0000705	4.82e-05	0.0002853	0.0127902	0.0302545	-0.0419399	0.0911629
RVE	Fisher	0.25	-0.0089165	7.30e-06	-0.0121348	-0.0056306	-0.0356660	0.0108344	-0.0485391	-0.0225224	0.0116499	0.0030605	0.0069931	0.0170182	0.0001449	0.0000720	4.89e-05	0.0002896	0.0132855	0.0303535	-0.0416434	0.0917111
RVE	Fisher TF	0.25	-0.0089388	9.20e-06	-0.0134198	-0.0044227	-0.0357552	0.0121479	-0.0536793	-0.0176909	0.0116504	0.0031698	0.0057559	0.0177705	0.0001456	0.0000743	3.31e-05	0.0003158	0.0003851	0.0253995	-0.0723426	0.0523955
RVE	Hedges	0.25	-0.0086741	6.90e-06	-0.0117676	-0.0054521	-0.0346965	0.0105305	-0.0470704	-0.0218084	0.0114562	0.0030115	0.0068941	0.0167275	0.0001401	0.0000698	4.75e-05	0.0002798	0.0122270	0.0302147	-0.0425272	0.0901897
RVE	Smith	0.25	-0.0091580	7.80e-06	-0.0128584	-0.0056401	-0.0366319	0.0111526	-0.0514337	-0.0225603	0.0118476	0.0031086	0.0072980	0.0173469	0.0001498	0.0000743	5.33e-05	0.0003009	0.0143791	0.0307191	-0.0382305	0.0957725
RVE	Swiger	0.25	-0.0089112	7.30e-06	-0.0121212	-0.0056289	-0.0356449	0.0108223	-0.0484848	-0.0225155	0.0116456	0.0030585	0.0069917	0.0170112	0.0001448	0.0000719	4.89e-05	0.0002894	0.0132648	0.0303513	-0.0416562	0.0916901
REML	Donner	0.50	-0.0046887	2.60e-06	-0.0070313	-0.0028291	-0.0093774	0.0032041	-0.0140626	-0.0056582	0.0101298	0.0032817	0.0055367	0.0166666	0.0001132	0.0000722	3.07e-05	0.0002778	-0.0225623	0.0203117	-0.0637735	0.0272628
REML	Fisher	0.50	-0.0047143	2.60e-06	-0.0070271	-0.0028216	-0.0094286	0.0032176	-0.0140541	-0.0056432	0.0101465	0.0032853	0.0055342	0.0166760	0.0001135	0.0000724	3.06e-05	0.0002781	-0.0230695	0.0203648	-0.0647129	0.0267452
REML	Fisher TF	0.50	-0.0119975	2.09e-05	-0.0183252	-0.0056447	-0.0239950	0.0091452	-0.0366503	-0.0112893	0.0152137	0.0044832	0.0076579	0.0232952	0.0002511	0.0001361	5.86e-05	0.0005427	-0.0005446	0.0220473	-0.0385554	0.0643555
REML	Hedges	0.50	-0.0046951	2.60e-06	-0.0070255	-0.0028390	-0.0093901	0.0031966	-0.0140509	-0.0056780	0.0101283	0.0032786	0.0055423	0.0166487	0.0001131	0.0000721	3.07e-05	0.0002772	-0.0221176	0.0202850	-0.0622546	0.0274322
REML	Smith	0.50	-0.0047318	2.60e-06	-0.0070290	-0.0028227	-0.0094637	0.0032358	-0.0140579	-0.0056454	0.0101635	0.0032911	0.0055352	0.0166775	0.0001139	0.0000726	3.06e-05	0.0002781	-0.0237785	0.0204854	-0.0647166	0.0266716
REML	Swiger	0.50	-0.0047139	2.60e-06	-0.0070271	-0.0028220	-0.0094278	0.0032170	-0.0140542	-0.0056440	0.0101461	0.0032851	0.0055343	0.0166753	0.0001135	0.0000724	3.06e-05	0.0002781	-0.0230428	0.0203606	-0.0646352	0.0267591
RVE	Donner	0.50	-0.0046818	2.60e-06	-0.0070300	-0.0028209	-0.0093635	0.0031993	-0.0140600	-0.0056417	0.0101308	0.0032846	0.0055340	0.0166750	0.0001132	0.0000723	3.06e-05	0.0002781	-0.0074028	0.0190845	-0.0471985	0.0403454
RVE	Fisher	0.50	-0.0047077	2.60e-06	-0.0070259	-0.0028127	-0.0094154	0.0032134	-0.0140518	-0.0056254	0.0101480	0.0032885	0.0055313	0.0166846	0.0001136	0.0000725	3.06e-05	0.0002784	-0.0073879	0.0190841	-0.0475849	0.0404717
RVE	Fisher TF	0.50	-0.0119960	2.09e-05	-0.0183229	-0.0056450	-0.0239920	0.0091421	-0.0366457	-0.0112900	0.0152124	0.0044821	0.0076581	0.0232928	0.0002511	0.0001360	5.86e-05	0.0005426	-0.0012317	0.0220094	-0.0396525	0.0637564
RVE	Hedges	0.50	-0.0046883	2.50e-06	-0.0070241	-0.0028373	-0.0093766	0.0031917	-0.0140481	-0.0056747	0.0101293	0.0032814	0.0055399	0.0166566	0.0001131	0.0000722	3.07e-05	0.0002774	-0.0073436	0.0190830	-0.0461815	0.0404359
RVE	Smith	0.50	-0.0047256	2.60e-06	-0.0070279	-0.0028138	-0.0094511	0.0032325	-0.0140558	-0.0056276	0.0101656	0.0032946	0.0055323	0.0166861	0.0001140	0.0000727	3.06e-05	0.0002784	-0.0073035	0.0191023	-0.0475624	0.0404760
RVE	Swiger	0.50	-0.0047073	2.60e-06	-0.0070260	-0.0028132	-0.0094145	0.0032127	-0.0140519	-0.0056264	0.0101476	0.0032883	0.0055314	0.0166839	0.0001136	0.0000725	3.06e-05	0.0002784	-0.0073856	0.0190835	-0.0475400	0.0404697
REML	Donner	0.90	0.0026873	8.00e-07	0.0015520	0.0042820	0.0029859	0.0009654	0.0017245	0.0047578	0.0049811	0.0014917	0.0026439	0.0083876	0.0000270	0.0000161	7.00e-06	0.0000704	-0.0482666	0.0255422	-0.1002644	0.0095753
REML	Fisher	0.90	0.0027883	9.00e-07	0.0015924	0.0044530	0.0030982	0.0010255	0.0017693	0.0049478	0.0050419	0.0015024	0.0026784	0.0084627	0.0000276	0.0000163	7.20e-06	0.0000716	-0.0489632	0.0256952	-0.1015855	0.0088566
REML	Fisher TF	0.90	-0.0029852	4.00e-06	-0.0064162	0.0003189	-0.0033169	0.0022094	-0.0071291	0.0003544	0.0054809	0.0017145	0.0024727	0.0088977	0.0000329	0.0000193	6.10e-06	0.0000792	-0.0077060	0.0240620	-0.0603324	0.0479065
REML	Hedges	0.90	0.0025850	7.00e-07	0.0015133	0.0040539	0.0028723	0.0009534	0.0016814	0.0045043	0.0049245	0.0014962	0.0026133	0.0082892	0.0000264	0.0000160	6.80e-06	0.0000687	-0.0477363	0.0256486	-0.0983524	0.0109008
REML	Smith	0.90	0.0029853	1.00e-06	0.0017077	0.0050108	0.0033170	0.0011057	0.0018974	0.0055676	0.0051620	0.0015152	0.0026858	0.0087250	0.0000289	0.0000168	7.20e-06	0.0000761	-0.0501228	0.0258765	-0.1054508	0.0077634
REML	Swiger	0.90	0.0027842	8.00e-07	0.0015912	0.0044430	0.0030935	0.0010227	0.0017680	0.0049367	0.0050394	0.0015019	0.0026775	0.0084583	0.0000276	0.0000163	7.20e-06	0.0000715	-0.0489314	0.0256892	-0.1015095	0.0088963
RVE	Donner	0.90	0.0026360	7.00e-07	0.0015095	0.0041099	0.0029289	0.0009001	0.0016772	0.0045665	0.0049578	0.0014872	0.0026054	0.0083483	0.0000267	0.0000160	6.80e-06	0.0000697	-0.0621029	0.0269319	-0.1215069	0.0034897
RVE	Fisher	0.90	0.0027390	7.00e-07	0.0015517	0.0042677	0.0030433	0.0009601	0.0017241	0.0047419	0.0050183	0.0014925	0.0026413	0.0084130	0.0000274	0.0000162	7.00e-06	0.0000708	-0.0612664	0.0273264	-0.1211787	0.0040406
RVE	Fisher TF	0.90	-0.0029847	4.00e-06	-0.0064147	0.0003188	-0.0033164	0.0022088	-0.0071275	0.0003542	0.0054806	0.0017143	0.0024727	0.0088978	0.0000329	0.0000193	6.10e-06	0.0000792	-0.0082505	0.0240534	-0.0601349	0.0468299
RVE	Hedges	0.90	0.0025347	6.00e-07	0.0014704	0.0039010	0.0028163	0.0008872	0.0016338	0.0043345	0.0049028	0.0014949	0.0025746	0.0082672	0.0000262	0.0000159	6.60e-06	0.0000683	-0.0626932	0.0265652	-0.1216705	0.0031726
RVE	Smith	0.90	0.0029397	9.00e-07	0.0017278	0.0047894	0.0032663	0.0010429	0.0019198	0.0053215	0.0051383	0.0014967	0.0026489	0.0086381	0.0000286	0.0000165	7.00e-06	0.0000746	-0.0596093	0.0284146	-0.1197092	0.0041511
RVE	Swiger	0.90	0.0027347	7.00e-07	0.0015506	0.0042585	0.0030386	0.0009573	0.0017229	0.0047317	0.0050157	0.0014922	0.0026403	0.0084092	0.0000273	0.0000162	7.00e-06	0.0000707	-0.0613019	0.0273074	-0.1212003	0.0040281

2 Pooled Estimate

2.1 Interactions

2.1.1 True ICC vs num pooled ICCs

perf_crit  |>  
  ggplot(aes(x = variance, y = rel_bias_icc_median, fill = model)) +
  geom_hline(yintercept = .05, linetype = "dashed")  + 
  geom_hline(yintercept = -.05, linetype = "dashed")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(var_combo_graph_f ~ icc_est_n, 
             labeller = label_bquote(rows = rho ==.(as.character(var_combo_graph_f)),
                                     cols = k ==.(icc_est_n))) + 
  labs(x = NULL, y = "RPB", fill = NULL) + 
  theme_bw() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

ggsave("RPB_numICC.png", width = 6.8, height = 5, units = "in")

perf_crit  |>  
  ggplot(aes(x = variance, y = rel_bias_icc_alt, fill = model)) + 
   geom_hline(yintercept = 1, linetype = "solid")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(var_combo_graph_f ~ icc_est_n, 
             labeller = label_bquote(rows = rho ==.(as.character(var_combo_graph_f)), 
                                            cols = k ==.(icc_est_n))) + 
  labs(x = NULL, y = "RPB - ratio", fill = NULL) + 
  theme_bw() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

ggsave("RPB_ratio_numICC.png", width = 6.8, height = 5, units = "in")

# perf_crit_MCSE <- perf_crit |> 
#   group_by(var_combo_graph, icc_est_n) |>  
#   summarise(mcse = max(bias_icc_mcse))


perf_crit  |>  
  ggplot(aes(x = variance, y = bias_icc_median, fill = model)) + 
   geom_hline(yintercept = 0, linetype = "solid", colour = "black")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(var_combo_graph_f ~ icc_est_n, 
             labeller = label_bquote(rows = rho ==.(as.character(var_combo_graph_f)), 
                                     cols = k ==.(icc_est_n))) + 
  labs(x = NULL, y = "PB", fill = NULL) + 
  theme_bw() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

ggsave("PB_numICC.png", width = 6.8, height = 5, units = "in")


# ggplot(perf_crit, aes(var_combo_graph_f, bias_icc, col=model ) ) +
#   facet_grid(variance ~ icc_est_n, 
#              labeller = label_bquote(rows = .(variance), 
#                                             cols = k ==.(icc_est_n))) +
#   geom_point( alpha=0.10,
#               position = position_dodge(width = 0.04) ) +
#   geom_smooth( se=FALSE ) + 
#   geom_hline( yintercept = 0 ) +
#   theme_minimal()


# ggplot(perf_crit, aes(var_combo_graph_f, bias_icc, col=skewness ) ) +
#   facet_grid(variance ~ icc_est_n, 
#              labeller = label_bquote(rows = .(variance), 
#                                             cols = k ==.(icc_est_n))) +
#   geom_point( alpha=0.10,
#               position = position_dodge(width = 0.04) ) +
#   geom_smooth( se=FALSE ) + 
#   geom_hline( yintercept = 0 ) +
#   theme_minimal()


perf_crit  |>  
  ggplot(aes(x = variance, y = RMSE, fill = model)) + 
  geom_hline(yintercept = .00, linetype = "dashed")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(var_combo_graph_f ~ icc_est_n, 
             labeller = label_bquote(rows = rho ==.(as.character(var_combo_graph_f)), 
                                            cols = k ==.(icc_est_n))) + 
  labs(x = NULL, y = "RMSE", fill = NULL) + 
  theme_bw() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

ggsave("RMSE_numICC.png", width = 6.8, height = 5, units = "in")

2.1.2 True ICC vs tau (between study heterogeneity)

perf_crit  |> 
  ggplot(aes(x = variance, y = rel_bias_icc, fill = model)) + 
   geom_hline(yintercept = .05, linetype = "dashed")  + 
   geom_hline(yintercept = -.05, linetype = "dashed")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(var_combo_graph_f ~ tau, 
             labeller = label_bquote(rows = rho ==.(as.character(var_combo_graph_f)), 
                                            cols = tau==.(tau))) + 
  labs(x = NULL, y = "RPB", fill = NULL) + 
  theme_bw() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

ggsave("RPB_tau.png", width = 6.8, height = 5, units = "in")

perf_crit  |> 
  ggplot(aes(x = variance, y = bias_icc, fill = model)) + 
   geom_hline(yintercept = 0, linetype = "solid")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(var_combo_graph_f ~ tau, 
             labeller = label_bquote(rows = rho ==.(as.character(var_combo_graph_f)), 
                                            cols = tau==.(tau))) + 
  labs(x = NULL, y = "PB", fill = NULL) + 
  theme_bw() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

ggsave("PB_tau.png", width = 6.8, height = 5, units = "in")



perf_crit  |> 
  ggplot(aes(x = variance, y = RMSE, fill = model)) + 
  geom_hline(yintercept = .00, linetype = "dashed") + 
  geom_boxplot(alpha = .5) + 
  facet_grid(var_combo_graph_f ~ tau, 
             labeller = label_bquote(rows = rho ==.(as.character(var_combo_graph_f)), 
                                            cols = tau==.(tau))) + 
  labs(x = NULL, y = "RMSE", fill = NULL) + 
  theme_bw() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

ggsave("RMSE_tau.png", width = 6.8, height = 5, units = "in")

2.1.3 True ICC vs number of clusters

perf_crit  |> 
  ggplot(aes(x = variance, y = rel_bias_icc, fill = model)) + 
   geom_hline(yintercept = .05, linetype = "dashed")  + 
   geom_hline(yintercept = -.05, linetype = "dashed")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(var_combo_graph_f ~ nj_size, 
             labeller = label_bquote(rows = rho ==.(as.character(var_combo_graph_f)), 
                                     cols = j== .(nj_size))) + 
  labs(x = NULL, y = "RPB", fill = NULL) + 
  theme_bw() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

ggsave("RPB_j.png", width = 6.8, height = 5, units = "in")

perf_crit  |> 
  ggplot(aes(x = variance, y = bias_icc, fill = model)) + 
   geom_hline(yintercept = 0, linetype = "solid")  + 
  # geom_hline(yintercept = -.05, linetype = "dashed")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid( var_combo_graph_f ~ nj_size, 
             labeller = label_bquote(rows = rho ==.(as.character(var_combo_graph_f)), 
                                     cols = j== .(nj_size))) + 
  labs(x = NULL, y = "PB", fill = NULL) + 
  theme_bw() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

ggsave("PB_j.png", width = 6.8, height = 5, units = "in")


perf_crit  |> 
  ggplot(aes(x = variance, y = RMSE, fill = model)) + 
  geom_hline(yintercept = .00, linetype = "dashed")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(var_combo_graph_f ~ nj_size, 
             labeller = label_bquote(rows = rho ==.(as.character(var_combo_graph_f)), 
                                     cols = j== .(nj_size))) + 
  labs(x = NULL, y = "RMSE", fill = NULL) + 
  theme_bw() +
  theme(
    legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

ggsave("RMSE_j.png", width = 6.8, height = 5, units = "in")

2.1.4 True ICC vs cluster size

perf_crit  |> 
  ggplot(aes(x = variance, y = rel_bias_icc, fill = model)) + 
   geom_hline(yintercept = .05, linetype = "solid")  + 
   geom_hline(yintercept = -.05, linetype = "dashed")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(var_combo_graph_f ~ n_bar_size, 
             labeller = label_bquote(rows = rho ==.(as.character(var_combo_graph_f)),
                                     cols = bar(n_j)==.(n_bar_size))) + 
  labs(x = NULL, y = "RPB", fill = NULL) + 
  theme_bw() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

ggsave("RPB_nj.png", width = 6.8, height = 5, units = "in")

perf_crit  |> 
  ggplot(aes(x = variance, y = bias_icc, fill = model)) + 
   geom_hline(yintercept = 0, linetype = "solid")  + 
  # geom_hline(yintercept = -.05, linetype = "dashed")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(var_combo_graph_f ~ n_bar_size, 
             labeller = label_bquote(rows = rho ==.(as.character(var_combo_graph_f)),
                                     cols = bar(n_j)==.(n_bar_size))) + 
  labs(x = NULL, y = "PB", fill = NULL) + 
  theme_bw() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

ggsave("PB_nj.png", width = 6.8, height = 5, units = "in")


perf_crit  |> 
  ggplot(aes(x = variance, y = RMSE, fill = model)) + 
  geom_hline(yintercept = .00, linetype = "dashed")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(var_combo_graph_f ~ n_bar_size, 
             labeller = label_bquote(rows = rho ==.(as.character(var_combo_graph_f)),
                                     cols = bar(n_j)==.(n_bar_size))) + 
  labs(x = NULL, y = "RMSE", fill = NULL) + 
  theme_bw() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

ggsave("RMSE_nj.png", width = 6.8, height = 5, units = "in")

2.1.5 True ICC vs Degree of unbalance

perf_crit  |> 
  ggplot(aes(x = variance, y = rel_bias_icc, fill = model)) + 
   geom_hline(yintercept = .05, linetype = "dashed")  + 
   geom_hline(yintercept = -.05, linetype = "dashed")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(var_combo_graph_f ~ n_bar_prop, 
             labeller = label_bquote(rows = rho ==.(as.character(var_combo_graph_f)),
                                     cols = zeta ==.(n_bar_prop))) + 
  labs(x = NULL, y = "RPB", fill = NULL) + 
  theme_bw() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

#ggsave("RPB_n_bar_prop.png", width = 6.8, height = 5, units = "in")

perf_crit  |> 
  ggplot(aes(x = variance, y = bias_icc, fill = model)) + 
  geom_hline(yintercept = 0, linetype = "solid")  + 
  # geom_hline(yintercept = .05, linetype = "dashed")  + 
  # geom_hline(yintercept = -.05, linetype = "dashed")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid( var_combo_graph_f ~ n_bar_prop, 
             labeller = label_bquote(rows = rho ==.(as.character(var_combo_graph_f)),
                                     cols = zeta ==.(n_bar_prop))) + 
  labs(x = NULL, y = "PB", fill = NULL) + 
  theme_bw() +
  theme(
    legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

#ggsave("PB_n_bar_prop.png", width = 6.8, height = 5, units = "in")


perf_crit  |> 
  ggplot(aes(x = variance, y = RMSE, fill = model)) + 
  geom_hline(yintercept = .00, linetype = "dashed")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid( var_combo_graph_f ~ n_bar_prop, 
             labeller = label_bquote(rows = rho ==.(as.character(var_combo_graph_f)),
                                     cols =zeta ==.(n_bar_prop))) + 
  labs(x = NULL, y = "RMSE", fill = NULL) + 
  theme_bw() +
  theme(
    legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

#ggsave("RMSE_n_bar_prop.png", width = 6.8, height = 5, units = "in")

2.1.6 number clusters vs cluster size

This is how the variance formulae perform across all generating ICC values.

perf_crit  |> 
  ggplot(aes(x = variance, y = rel_bias_icc, fill = model)) + 
  geom_hline(yintercept = .05, linetype = "dashed")  + 
  geom_hline(yintercept = -.05, linetype = "dashed")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(nj_size ~ n_bar_size, 
             labeller = label_bquote(rows = j ==.(nj_size), 
                                     cols = bar(n_j)==.(n_bar_size))) + 
  labs(x = NULL, y = "RPB", fill = NULL) + 
  theme_bw() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

#ggsave("RPB_j_nj.png", width = 6.8, height = 5, units = "in")

perf_crit  |> 
  ggplot(aes(x = variance, y = bias_icc, fill = model)) + 
  geom_hline(yintercept = 0, linetype = "solid")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(nj_size ~ n_bar_size, 
             labeller = label_bquote(rows = j ==.(nj_size), 
                                     cols = bar(n_j)==.(n_bar_size))) + 
  labs(x = NULL, y = "PB", fill = NULL) + 
  theme_bw() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

#ggsave("PB_j_nj.png", width = 6.8, height = 5, units = "in")

perf_crit  |> ggplot(aes(x = variance, y = RMSE, fill = model)) + 
  geom_hline(yintercept = .00, linetype = "dashed")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(nj_size ~ n_bar_size, 
             labeller = label_bquote(rows = j ==.(nj_size), 
                                            cols = bar(n_j)==.(n_bar_size))) + 
  labs(x = NULL, y = "RMSE", fill = NULL) + 
  theme_bw() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

#ggsave("RMSE_j_nj.png", width = 6.8, height = 5, units = "in")

perf_crit  |> 
  ggplot(aes(x = variance, y = rel_bias_icc, fill = nj_size)) + 
  geom_hline(yintercept = .05, linetype = "dashed")  + 
  geom_hline(yintercept = -.05, linetype = "dashed")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(var_combo_graph_f ~ n_bar_size, 
             labeller = label_bquote(rows = rho ==.(as.character(var_combo_graph_f)), 
                                     cols = bar(n_j)==.(n_bar_size))) + 
  labs(x = NULL, y = "RPB", fill = NULL) + 
  theme_bw() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

ggsave("RPB_j_nj_nomodel.png", width = 6.8, height = 5, units = "in")

perf_crit  |> 
  ggplot(aes(x = variance, y = bias_icc, fill = nj_size)) + 
  geom_hline(yintercept = 0, linetype = "solid")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(var_combo_graph_f ~ n_bar_size, 
             labeller = label_bquote(rows = rho ==.(as.character(var_combo_graph_f)), 
                                     cols = bar(n_j)==.(n_bar_size))) + 
  labs(x = NULL, y = "PB", fill = NULL) + 
  theme_bw() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

ggsave("PB_j_nj_nomodel.png", width = 6.8, height = 5, units = "in")

perf_crit  |> ggplot(aes(x = variance, y = RMSE, fill = nj_size)) + 
  geom_hline(yintercept = .00, linetype = "dashed")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(var_combo_graph_f ~ n_bar_size, 
             labeller = label_bquote(rows = rho ==.(as.character(var_combo_graph_f)),
                                     cols = bar(n_j)==.(n_bar_size))) + 
  labs(x = NULL, y = "RMSE", fill = NULL) + 
  theme_bw() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

#ggsave("RMSE_j_nj.png", width = 6.8, height = 5, units = "in")

perf_crit  |> 
  ggplot(aes(x = variance, y = rel_bias_icc, fill = model)) + 
  geom_hline(yintercept = .05, linetype = "dashed")  + 
  geom_hline(yintercept = -.05, linetype = "dashed")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(n_bar_prop ~ n_bar_size, 
             labeller = label_bquote(rows = zeta ==.(n_bar_prop), 
                                            cols = bar(n_j)==.(n_bar_size))) + 
  labs(x = NULL, y = "RPB", fill = NULL) + 
  theme_bw() +
  theme(
    legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),

        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

perf_crit  |> 
  ggplot(aes(x = variance, y = bias_icc, fill = model)) + 
  geom_hline(yintercept = 0, linetype = "solid")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(n_bar_prop ~ n_bar_size, 
             labeller = label_bquote(rows = zeta ==.(n_bar_prop), 
                                            cols = bar(n_j)==.(n_bar_size))) + 
  labs(x = NULL, y = "PB", fill = NULL) + 
  theme_bw() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

filtered \(\rho = 0.50\)

perf_crit  |> filter(var_combo_graph == 0.50) |> 
  ggplot(aes(x = variance, y = rel_bias_icc, fill = model)) + 
  geom_hline(yintercept = .05, linetype = "dashed")  + 
  geom_hline(yintercept = -.05, linetype = "dashed")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(nj_size ~ n_bar_size, 
             labeller = label_bquote(rows = j ==.(nj_size),
                                            cols = bar(n_j)==.(n_bar_size))) + 
  labs(x = NULL, y = "RPB", fill = NULL) + 
  theme_bw() +
  theme(
    legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),

        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

perf_crit  |> filter(var_combo_graph == 0.50) |> 
  ggplot(aes(x = variance, y = bias_icc, fill = model)) + 
  geom_hline(yintercept = 0, linetype = "solid")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(nj_size ~ n_bar_size, 
             labeller = label_bquote(rows = j ==.(nj_size),
                                            cols = bar(n_j)==.(n_bar_size))) + 
  labs(x = NULL, y = "PB", fill = NULL) + 
  theme_bw() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

perf_crit  |> filter(var_combo_graph == 0.50) |> 
  ggplot(aes(x = variance, y = rel_bias_icc, fill = model)) + 
  geom_hline(yintercept = .05, linetype = "dashed")  + 
  geom_hline(yintercept = -.05, linetype = "dashed")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(n_bar_prop ~ n_bar_size, 
             labeller = label_bquote(rows = zeta ==.(n_bar_prop),
                                            cols = bar(n_j)==.(n_bar_size))) + 
  labs(x = NULL, y = "RPB", fill = NULL) + 
  theme_bw() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

perf_crit  |> filter(var_combo_graph == 0.50) |> 
  ggplot(aes(x = variance, y = bias_icc, fill = model)) + 
  geom_hline(yintercept = 0, linetype = "solid")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(n_bar_prop ~ n_bar_size, 
             labeller = label_bquote(rows = zeta ==.(n_bar_prop), 
                                            cols = bar(n_j)==.(n_bar_size))) + 
  labs(x = NULL, y = "PB", fill = NULL) + 
  theme_bw() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

filtered \(\rho = 0.90\)

perf_crit  |> filter(var_combo_graph == 0.90) |> 
  ggplot(aes(x = variance, y = rel_bias_icc, fill = model)) + 
  geom_hline(yintercept = .05, linetype = "dashed")  + 
  geom_hline(yintercept = -.05, linetype = "dashed")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(nj_size ~ n_bar_size, 
             labeller = label_bquote(rows = j ==.(nj_size),
                                            cols = bar(n_j)==.(n_bar_size))) + 
  labs(x = NULL, y = "RPB", fill = NULL) + 
  theme_bw() +
  theme(
    legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),

        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

perf_crit  |> filter(var_combo_graph == 0.90) |> 
  ggplot(aes(x = variance, y = bias_icc, fill = model)) + 
  geom_hline(yintercept = 0, linetype = "solid")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(nj_size ~ n_bar_size, 
             labeller = label_bquote(rows = j ==.(nj_size),
                                            cols = bar(n_j)==.(n_bar_size))) + 
  labs(x = NULL, y = "PB", fill = NULL) + 
  theme_bw() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

perf_crit  |> filter(var_combo_graph == 0.90) |> 
  ggplot(aes(x = variance, y = rel_bias_icc, fill = model)) + 
  geom_hline(yintercept = .05, linetype = "dashed")  + 
  geom_hline(yintercept = -.05, linetype = "dashed")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(n_bar_prop ~ n_bar_size, 
             labeller = label_bquote(rows = zeta ==.(n_bar_prop), 
                                            cols = bar(n_j)==.(n_bar_size))) + 
  labs(x = NULL, y = "RPB", fill = NULL) + 
  theme_bw() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

perf_crit  |> filter(var_combo_graph == 0.90) |> 
  ggplot(aes(x = variance, y = bias_icc, fill = model)) + 
  geom_hline(yintercept = 0, linetype = "solid")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(n_bar_prop ~ n_bar_size, 
             labeller = label_bquote(rows = zeta ==.(n_bar_prop), 
                                            cols = bar(n_j)==.(n_bar_size))) + 
  labs(x = NULL, y = "PB", fill = NULL) + 
  theme_bw() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

filtered \(\rho = 0.05\)

perf_crit  |> filter(var_combo_graph == 0.05) |> 
  ggplot(aes(x = variance, y = rel_bias_icc, fill = model)) + 
  geom_hline(yintercept = .05, linetype = "dashed")  + 
  geom_hline(yintercept = -.05, linetype = "dashed")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(nj_size ~ n_bar_size, 
             labeller = label_bquote(rows = j ==.(nj_size),
                                            cols = bar(n_j)==.(n_bar_size))) + 
  labs(x = NULL, y = "RPB", fill = NULL) + 
  theme_bw() +
  theme(
    legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),

        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

perf_crit  |> filter(var_combo_graph == 0.05) |> 
  ggplot(aes(x = variance, y = bias_icc, fill = model)) + 
  geom_hline(yintercept = 0, linetype = "solid")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(nj_size ~ n_bar_size, 
             labeller = label_bquote(rows = j ==.(nj_size),
                                            cols = bar(n_j)==.(n_bar_size))) + 
  labs(x = NULL, y = "PB", fill = NULL) + 
  theme_bw() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

perf_crit  |> filter(var_combo_graph == 0.05) |> 
  ggplot(aes(x = variance, y = rel_bias_icc, fill = nj_size)) + 
  geom_hline(yintercept = .05, linetype = "dashed")  + 
  geom_hline(yintercept = -.05, linetype = "dashed")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(n_bar_prop ~ n_bar_size, 
             labeller = label_bquote(rows = zeta ==.(n_bar_prop), 
                                            cols = bar(n_j)==.(n_bar_size))) + 
  labs(x = NULL, y = "RPB", fill = NULL) + 
  theme_bw() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

perf_crit  |> filter(var_combo_graph == 0.05) |> 
  ggplot(aes(x = variance, y = bias_icc, fill = nj_size)) + 
  geom_hline(yintercept = 0, linetype = "solid")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(n_bar_prop ~ n_bar_size, 
             labeller = label_bquote(rows = zeta ==.(n_bar_prop), 
                                            cols = bar(n_j)==.(n_bar_size))) + 
  labs(x = NULL, y = "PB", fill = NULL) + 
  theme_bw() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

3 SE of ICC

3.1 Interactions

3.1.1 True ICC vs num pooled ICCs

perf_crit  |> ggplot(aes(x = variance, y = rel_bias_se_ICC, fill = model)) + 
  geom_hline(yintercept = .05, linetype = "dashed")  + 
  geom_hline(yintercept = -.05, linetype = "dashed")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(var_combo_graph_f ~ icc_est_n, 
             labeller = label_bquote(rows = rho ==.(as.character(var_combo_graph_f)),
                                     cols = k ==.(icc_est_n))) + 
  labs(x = NULL, y = "RSEB", fill = NULL) + 
  theme_bw() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

ggsave("RSEB_numICC.png", width = 6.8, height = 5, units = "in")

# perf_crit  |> ggplot(aes(x = variance, y = rel_bias_se_ICC, fill = model)) + 
#   geom_hline(yintercept = .05, linetype = "dashed")  + 
#   geom_hline(yintercept = -.05, linetype = "dashed")  + 
#   geom_boxplot(alpha = .5) + 
#   facet_grid( var_combo_graph_f ~ icc_est_n, 
#              labeller = label_bquote(rows = rho ==.(as.character(var_combo_graph_f)), 
#                                             cols = k ==.(icc_est_n))
#              ) + 
#   labs(x = NULL, y = "RSEB", fill = NULL) + 
#   theme_bw() +
#   theme(
#     legend.position = "bottom",
#         axis.text.x = element_text(size = 6.5),
#         axis.text.y = element_text(size = 6.5),
# 
#         plot.caption=element_text(hjust = 0, size = 6),
#        # legend.title=element_text(size = 6), 
#         legend.text=element_text(size = 6.5),
#         axis.title.x = element_text(size=6.5),
#         axis.title.y = element_text(size=6.5),
#       #  strip.background =element_rect(fill="white")
#         strip.text.x = element_text(size = 6.5, face="bold"),
#         strip.text.y = element_text(size = 6.5, face="bold")
#       )

#ggsave("RPB_se_numICC.png", width = 6.8, height = 5, units = "in")

perf_crit |>
  ggplot(aes(x = as.character(variance), y = rel_bias_se_ICC, 
             color = model, 
             shape = as.character(tau)
             )
         ) +
  geom_point(alpha = .8, 
           #  position = position_jitter(height = .3), 
             size = 1)   + 
  geom_hline(yintercept = .05, linetype = "dashed")  + 
  geom_hline(yintercept = -.05, linetype = "dashed")  +
  facet_grid( var_combo_graph_f ~ icc_est_n, 
             labeller = label_bquote(rows = rho ==.(as.character(var_combo_graph_f)), 
                                            cols = k ==.(icc_est_n))) + 
  labs(
    x = NULL, 
    y ="RSEB",
    color = "", shape = ""
  ) + 
  theme_bw() +
  theme(
    plot.caption=element_text(hjust = 0, size = 10),
    legend.position= "bottom",
    panel.spacing.x = unit(5, "mm"), 
    panel.spacing.y = unit(5, "mm"),
  )

ggsave("RSEB_numICC_alt2.png", width = 6.8, height = 5, units = "in")

3.1.2 True ICC vs tau (between study )

perf_crit  |> ggplot(aes(x = variance, y = rel_bias_se_ICC, fill = model)) + 
  geom_hline(yintercept = .05, linetype = "dashed")  + 
  geom_hline(yintercept = -.05, linetype = "dashed")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid( var_combo_graph_f ~ tau, 
             labeller = label_bquote(rows = rho ==.(as.character(var_combo_graph_f)), 
                                            cols = tau==.(tau))) + 
  labs(x = NULL, y = "RSEB", fill = NULL) + 
  theme_bw() +
  theme(
    legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),

        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
      #  strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

ggsave("RSEB_tau.png", width = 6.8, height = 5, units = "in")

3.1.3 True ICC vs number of clusters

perf_crit  |> ggplot(aes(x = variance, y = rel_bias_se_ICC, fill = model)) + 
  geom_hline(yintercept = .05, linetype = "dashed")  + 
  geom_hline(yintercept = -.05, linetype = "dashed")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid( var_combo_graph_f ~ nj_size, 
             labeller = label_bquote(rows = rho ==.(as.character(var_combo_graph_f)), 
                                            cols = j== .(nj_size))) + 
  labs(x = NULL, y = "RSEB", fill = NULL) + 
  theme_bw() +
  theme(
    legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),

        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
      #  strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

#ggsave("RPB_se_j.png", width = 6.8, height = 5, units = "in")

3.1.4 True ICC vs cluster size

TASHA – interestingly for SE of the pooled estimate, despite the bias in the pooled estimate.

perf_crit  |> ggplot(aes(x = variance, y = rel_bias_se_ICC, fill = model)) + 
  geom_hline(yintercept = .05, linetype = "dashed")  + 
  geom_hline(yintercept = -.05, linetype = "dashed")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid( var_combo_graph_f ~ n_bar_size, 
             labeller = label_bquote(rows = rho ==.(as.character(var_combo_graph_f)), 
                                            cols = bar(n_j)==.(n_bar_size))) + 
  labs(x = NULL, y = "RSEB", fill = NULL) + 
  theme_bw() +
  theme(
    legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),

        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
      #  strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

ggsave("RSEB_nj.png", width = 6.8, height = 5, units = "in")

3.1.5 True ICC vs Degree of unbalance

perf_crit  |> 
  ggplot(aes(x = variance, y = rel_bias_se_ICC, fill = model)) + 
  geom_hline(yintercept = .05, linetype = "dashed")  + 
  geom_hline(yintercept = -.05, linetype = "dashed")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid( var_combo_graph_f ~ n_bar_prop, 
             labeller = label_bquote(rows = rho ==.(as.character(var_combo_graph_f)), 
                                            cols = zeta ==.(n_bar_prop))) + 
  labs(x = NULL, y = "RSEB", fill = NULL) + 
  theme_bw() +
  theme(
    legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),

        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
      #  strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

#ggsave("RSEB_n_bar_prop.png", width = 6.8, height = 5, units = "in")

3.1.6 number clusters vs cluster size

perf_crit  |> ggplot(aes(x = variance, y = rel_bias_se_ICC, fill = model)) + 
  geom_hline(yintercept = .05, linetype = "dashed")  + 
  geom_hline(yintercept = -.05, linetype = "dashed")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid( nj_size ~ n_bar_size, 
             labeller = label_bquote(rows = j ==.(nj_size), 
                                            cols = bar(n_j)==.(n_bar_size))) + 
  labs(x = NULL, y = "RSEB", fill = NULL) + 
  theme_bw() +
  theme(
    legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),

        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
      #  strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

#ggsave("RSEB_j_nj.png", width = 6.8, height = 5, units = "in")

4 L2 Variance

4.1 Main Effects

4.1.1 L2 variance

perf_crit2  |> ggplot(aes(x = factor(var_combo_graph), y = rel_bias_clus)) + 
  geom_hline(yintercept = .05, linetype = "dashed")  + 
  geom_hline(yintercept = -.05, linetype = "dashed")  + 
  geom_boxplot(alpha = .5) + 
 # facet_grid( rows = vars(var_combo_graph_f)
              #labeller = label_parsed
#              ) + 
  labs(x = "ICC Parameter", y = "RPB") + 
  theme_bw() +
  theme(legend.position = "none",
        axis.text.x = element_text(size = 6, face="bold"),
        plot.caption=element_text(hjust = 0, size = 8),
        legend.title=element_text(size = 8), 
        legend.text=element_text(size = 8),
        axis.title.x = element_text(size=8),
        axis.title.y = element_text(size=8),
        strip.background =element_rect(fill="white"))

perf_crit2  |> ggplot(aes(x = factor(var_combo_graph), y = rel_bias_clus_ratio)) + 
  geom_hline(yintercept = 1, linetype = "solid")  + 
  geom_boxplot(alpha = .5) + 
 # facet_grid( rows = vars(var_combo_graph_f)
              #labeller = label_parsed
#              ) + 
  labs(x = "ICC Parameter", y = "RPB Ratio") + 
  theme_bw() +
  theme(legend.position = "none",
        axis.text.x = element_text(size = 6, face="bold"),
        plot.caption=element_text(hjust = 0, size = 8),
        legend.title=element_text(size = 8), 
        legend.text=element_text(size = 8),
        axis.title.x = element_text(size=8),
        axis.title.y = element_text(size=8),
        strip.background =element_rect(fill="white"))

4.1.2 L2 variance by tau

perf_crit2  |> ggplot(aes(x = factor(var_combo_graph), y = rel_bias_clus)) + 
  geom_hline(yintercept = .05, linetype = "dashed")  + 
  geom_hline(yintercept = -.05, linetype = "dashed")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(cols = vars(tau )
             , 
             labeller = label_bquote(cols = tau ==.(tau))) + 
  labs(x = "ICC Parameter", y = "RPB") + 
  theme_bw() +
  theme(legend.position = "none",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
      #  plot.caption=element_text(hjust = 0, size = 8),
      #  legend.title=element_text(size = 8), 
      #  legend.text=element_text(size = 8),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

ggsave("RPB_L2_tau_mean.png", width = 6.8, height = 5, units = "in")



perf_crit2  |> ggplot(aes(x = factor(var_combo_graph), y = rel_bias_clus_ratio)) + 
  geom_hline(yintercept = 1, linetype = "solid")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(cols = vars(tau )
             , 
             labeller = label_bquote(cols = tau ==.(tau))) + 
  labs(x = "ICC Parameter", y = "RPB Ratio") + 
  theme_bw() +
  theme(legend.position = "none",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
      #  plot.caption=element_text(hjust = 0, size = 8),
      #  legend.title=element_text(size = 8), 
      #  legend.text=element_text(size = 8),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

4.1.3 L2 variance by number of clusters

perf_crit2  |> ggplot(aes(x = factor(var_combo_graph), y = rel_bias_clus)) + 
  geom_hline(yintercept = .05, linetype = "dashed")  + 
  geom_hline(yintercept = -.05, linetype = "dashed")  + 
#  geom_hline(yintercept = 0, linetype = "solid", colour = "red")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(cols = vars(nj_size )
             , 
             labeller = label_bquote(cols = j ==.(nj_size))) + 
  labs(x = "ICC Parameter", y = "RPB") + 
  theme_bw() +
  theme(legend.position = "none",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
      #  plot.caption=element_text(hjust = 0, size = 8),
      #  legend.title=element_text(size = 8), 
      #  legend.text=element_text(size = 8),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

ggsave("RPB_L2_j_mean.png", width = 6.8, height = 5, units = "in")



perf_crit2  |> ggplot(aes(x = factor(var_combo_graph), y = rel_bias_clus_ratio)) + 
  geom_hline(yintercept = 1, linetype = "solid")  + 
 
  geom_boxplot(alpha = .5) + 
  facet_grid(cols = vars(nj_size )
             , 
             labeller = label_bquote(cols = j ==.(nj_size))) + 
  labs(x = "ICC Parameter", y = "RPB Ratio") + 
  theme_bw() +
  theme(legend.position = "none",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
      #  plot.caption=element_text(hjust = 0, size = 8),
      #  legend.title=element_text(size = 8), 
      #  legend.text=element_text(size = 8),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

4.1.4 L1 variance by clustersize

perf_crit2  |> ggplot(aes(x = factor(var_combo_graph), y = rel_bias_res)) + 
  geom_hline(yintercept = .05, linetype = "dashed")  + 
  geom_hline(yintercept = -.05, linetype = "dashed")  + 
#  geom_hline(yintercept = 0, linetype = "solid", colour = "red")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(cols = vars(n_bar_size )
             , 
             labeller = label_bquote(cols = nj ==.(n_bar_size))) + 
  labs(x = "ICC Parameter", y = "RPB") + 
  theme_bw() +
  theme(legend.position = "none",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
      #  plot.caption=element_text(hjust = 0, size = 8),
      #  legend.title=element_text(size = 8), 
      #  legend.text=element_text(size = 8),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

ggsave("RPB_L2_nj_mean.png", width = 6.8, height = 5, units = "in")



perf_crit2  |> ggplot(aes(x = factor(var_combo_graph), y = rel_bias_res_ratio)) + 
   geom_hline(yintercept = 1, linetype = "solid")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(cols = vars(n_bar_size )
             , 
             labeller = label_bquote(cols = nj ==.(n_bar_size))) + 
  labs(x = "ICC Parameter", y = "RPB Ratio") + 
  theme_bw() +
  theme(legend.position = "none",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
      #  plot.caption=element_text(hjust = 0, size = 8),
      #  legend.title=element_text(size = 8), 
      #  legend.text=element_text(size = 8),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

5 Evaluating the Distribution ICC estimates

Note for the Fisher TF pooled estimates, they have already been back transformed for all of these analyses.

load("data/results_conv.RData")
load("data/results_pooled_est.RData")

6 Distributions of the ICC Estimates from primary studies

skewness <- function(x) {
  
  variance = var(x)
  K = length(x)
  skewness = (1/(K*sqrt(variance)^3)*sum((x-mean(x))^3))
  
  return(skewness)
  
}

kurtosis <- function(x) {
  
  variance = var(x)
  K = length(x)
  kurtosis = (1/(K*variance^2))*sum((x-mean(x))^4)
  
  return(kurtosis)
  
}

First, which set of conditions resulted in the lowest raw bias and relative bias?

For \(\rho = 0.5\), it is when \(k = 50\), \(n_j \sim U[30, 50]\), \(\overline{n_j} \sim U[10, 30]\), \(n_{prop} = 0.1\), and \(\tau = 0.2\). I will use that set for the distributions.

# perf_crit |>  group_by(ICC_true) |>  slice(which.min(bias_icc)) |> select(ICC_true, method, var_icc_name, icc_est_n, nj_size, n_bar_size, n_bar_prop, tau)
# 
# perf_crit |>  group_by(ICC_true) |>  slice(which.min(rel_bias_icc)) |> select(ICC_true, method, var_icc_name, icc_est_n, nj_size, n_bar_size, n_bar_prop, tau)


results_all_converged_dist_icc_est <- results_all_converged |> filter(icc_est_n == 50) |>  filter(nj_size == "U[30, 50]") |>  filter(n_bar_size == "U[10, 30]") |>  filter(n_bar_prop == 0.1) |>  filter(tau == 0.02) 


# results_all_converged_dist_icc_est <- results_all_converged |> filter(nj_size == "U[30, 50]") |> filter(n_bar_size == "U[10, 30]") |> 
#   filter(n_bar_prop == 0.1) |> filter(tau == 0.02 ) |> filter(icc_est_n == 50) 

6.1 Distribution of ICC estimates when \(\rho = 0.05\):

dist_0.05 <- results_all_converged_dist_icc_est |> ungroup() |>  filter(var_combo_graph == 0.05) |> select(clustervar, residvar) |> distinct()  |>
  mutate(icc_list = map2(clustervar, residvar, ~ map2(.x, .y, ~ .x / (.x + .y)))) |>  select(icc_list) |> unnest(icc_list) |>  unlist() 

dist_0.05  <- data.frame(icc = dist_0.05)

ggplot(dist_0.05 , aes(x = icc)) +
  geom_histogram(fill = "#92D895", color = "#e9ecef", alpha = 0.6, binwidth = 0.01) +
  geom_vline(xintercept = 0.05, linetype = "solid") + 
  geom_vline(xintercept = mean(dist_0.05$icc), linetype = "dashed") + 
  geom_vline(xintercept = median(dist_0.05$icc), linetype = "dotted") + 
  labs(title = "Histogram of ICC Values for rho = .05 and cluster size 'U[10, 30]'", x = "ICC Estimate", y = "Frequency") +
  theme_minimal()

This distribution has a mean of 0.0500555, variance of 0.0010352, median of 0.0464974, a skewness of 0.5933149, and kurtosis of 3.1094325.

6.2 Distribution of ICC estimates when \(\rho = 0.10\):

dist_0.10 <- results_all_converged_dist_icc_est |> ungroup() |>  filter(var_combo_graph == 0.10) |> select(clustervar, residvar) |> distinct()  |>
  mutate(icc_list = map2(clustervar, residvar, ~ map2(.x, .y, ~ .x / (.x + .y)))) |>  select(icc_list) |> unnest(icc_list) |>  unlist() 

dist_0.10  <- data.frame(icc = dist_0.10)

ggplot(dist_0.10 , aes(x = icc)) +
  geom_histogram(fill = "#92D895", color = "#e9ecef", alpha = 0.6, binwidth = 0.01) +
      geom_vline(xintercept = 0.1, linetype = "solid") + 
  geom_vline(xintercept = mean(dist_0.10$icc), linetype = "dashed") + 
  geom_vline(xintercept = median(dist_0.10$icc), linetype = "dotted") + 
  labs(title = "Histogram of ICC Values for rho = .10 and cluster size 'U[10, 30]'", x = "ICC Estimate", y = "Frequency") +
  theme_minimal()

This distribution has a mean of 0.0993899, variance of 0.0015262, median of 0.0964493, a skewness of 0.4208593, and kurtosis of 3.1433641.

6.3 Distribution of ICC estimates when \(\rho = 0.15\):

dist_0.15 <- results_all_converged_dist_icc_est |> ungroup() |>  filter(var_combo_graph == 0.15) |> select(clustervar, residvar) |> distinct()  |>
  mutate(icc_list = map2(clustervar, residvar, ~ map2(.x, .y, ~ .x / (.x + .y)))) |>  select(icc_list) |> unnest(icc_list) |>  unlist() 

dist_0.15  <- data.frame(icc = dist_0.15)

ggplot(dist_0.15 , aes(x = icc)) +
  geom_histogram(fill = "#92D895", color = "#e9ecef", alpha = 0.6, binwidth = 0.01) +
      geom_vline(xintercept = 0.15, linetype = "solid") + 
  geom_vline(xintercept = mean(dist_0.15$icc), linetype = "dashed") + 
  geom_vline(xintercept = median(dist_0.15$icc), linetype = "dotted") + 
  labs(title = "Histogram of ICC Values for rho = .15 and cluster size 'U[10, 30]'", x = "ICC Estimate", y = "Frequency") +
  theme_minimal()

This distribution has a mean of 0.1482727, variance of 0.0020559, median of 0.1457706, a skewness of 0.3036363, and kurtosis of 3.0933228.

6.4 Distribution of ICC estimates when \(\rho = 0.25\):

dist_0.25 <- results_all_converged_dist_icc_est |> ungroup() |>  filter(var_combo_graph == 0.25) |> select(clustervar, residvar) |> distinct()  |>
  mutate(icc_list = map2(clustervar, residvar, ~ map2(.x, .y, ~ .x / (.x + .y)))) |>  select(icc_list) |> unnest(icc_list) |>  unlist() 

dist_0.25  <- data.frame(icc = dist_0.25 )

ggplot(dist_0.25 , aes(x = icc)) +
  geom_histogram(fill = "#92D895", color = "#e9ecef", alpha = 0.6, binwidth = 0.01) +
  geom_vline(xintercept = 0.25, linetype = "solid") + 
  geom_vline(xintercept = mean(dist_0.25$icc), linetype = "dashed") + 
  geom_vline(xintercept = median(dist_0.25$icc), linetype = "dotted") + 
  labs(title = "Histogram of ICC Values for rho = .25 and cluster size 'U[10, 30]'", x = "ICC Estimate", y = "Frequency") +
  theme_minimal()

This distribution has a mean of 0.2471755, variance of 0.0031014, median of 0.245689, a skewness of 0.140802, and kurtosis of 2.9323949.

6.5 Distribution of ICC estimates when \(\rho = 0.50\):

dist_0.5 <- results_all_converged_dist_icc_est |> ungroup() |>  filter(var_combo_graph == 0.5) |> select(clustervar, residvar) |> distinct()  |>
  mutate(icc_list = map2(clustervar, residvar, ~ map2(.x, .y, ~ .x / (.x + .y)))) |>  select(icc_list) |> unnest(icc_list) |>  unlist() 

dist_0.5  <- data.frame(icc = dist_0.5 )

ggplot(dist_0.5 , aes(x = icc)) +
 geom_histogram(fill = "#92D895", color = "#e9ecef", alpha = 0.6, binwidth = 0.01) +
  geom_vline(xintercept = 0.50, linetype = "solid") + 
  geom_vline(xintercept = mean(dist_0.5$icc), linetype = "dashed") + 
  geom_vline(xintercept = median(dist_0.5$icc), linetype = "dotted") + 
  labs(title = "Histogram of ICC Values for rho = .50 and cluster size 'U[10, 30]'", x = "ICC Estimate", y = "Frequency") +
  theme_minimal()

This distribution has a mean of 0.4923617, variance of 0.0044361, median of 0.4947604, a skewness of -0.2104669, and kurtosis of 3.0122678.

compare this to a distribution that has the larger average cluster sizes condition:

results_all_converged_dist_icc_est2 <- results_all_converged |> filter(nj_size == "U[30, 50]") |> filter(n_bar_size == "U[30, 50]") |> 
  filter(n_bar_prop == 0.1) |> filter(tau == 0.02   ) |> filter(icc_est_n == 50) 

dist_0.5 <- results_all_converged_dist_icc_est2 |> ungroup() |>  filter(var_combo_graph == 0.50) |> select(clustervar, residvar) |> distinct()  |>
  mutate(icc_list = map2(clustervar, residvar, ~ map2(.x, .y, ~ .x / (.x + .y)))) |>  select(icc_list) |> unnest(icc_list) |>  unlist() 

dist_0.5  <- data.frame(icc = dist_0.5)

ggplot(dist_0.5 , aes(x = icc)) +
  geom_histogram(fill = "#92D895", color = "#e9ecef", alpha = 0.6, binwidth = 0.01) +
  geom_vline(xintercept = 0.50, linetype = "solid") + 
  geom_vline(xintercept = mean(dist_0.5$icc), linetype = "dashed") + 
  geom_vline(xintercept = median(dist_0.5$icc), linetype = "dotted") + 
  labs(title = "Histogram of ICC Values for rho = .50 and cluster size 'U[30, 50]'", x = "ICC Estimate", y = "Frequency") +
  theme_minimal()

This distribution has a mean of 0.4928036, variance of 0.0040805, median of 0.4950266, a skewness of -0.1914525, and kurtosis of 2.9920759.

6.6 Distribution of ICC estimates when \(\rho = 0.90\):

dist_0.9 <- results_all_converged_dist_icc_est |> ungroup() |>  filter(var_combo_graph == 0.9) |> select(clustervar, residvar) |> distinct()  |>
  mutate(icc_list = map2(clustervar, residvar, ~ map2(.x, .y, ~ .x / (.x + .y)))) |>  select(icc_list) |> unnest(icc_list) |>  unlist() 

dist_0.9  <- data.frame(icc = dist_0.9)

ggplot(dist_0.9 , aes(x = icc)) +
  geom_histogram(fill = "#92D895", color = "#e9ecef", alpha = 0.6, binwidth = 0.01) +
  geom_vline(xintercept = 0.9, linetype = "solid") + 
  geom_vline(xintercept = mean(dist_0.9$icc), linetype = "dashed") + 
  geom_vline(xintercept = median(dist_0.9$icc), linetype = "dotted") + 
  labs(title = "Histogram of ICC Values for rho = .9 and cluster size 'U[10, 30]'", x = "ICC Estimate", y = "Frequency") +
  theme_minimal()

This distribution has a mean of 0.8956405, variance of 0.0011112, median of 0.8993834, a skewness of -0.691453, and kurtosis of 3.6654658.

dist_0.9 <- results_all_converged_dist_icc_est2 |> ungroup() |>  filter(var_combo_graph == 0.9) |> select(clustervar, residvar) |> distinct()  |>
  mutate(icc_list = map2(clustervar, residvar, ~ map2(.x, .y, ~ .x / (.x + .y)))) |>  select(icc_list) |> unnest(icc_list) |>  unlist() 

dist_0.9  <- data.frame(icc = dist_0.9)

ggplot(dist_0.9 , aes(x = icc)) +
  geom_histogram(fill = "#92D895", color = "#e9ecef", alpha = 0.6, binwidth = 0.01) +
  geom_vline(xintercept = 0.9, linetype = "solid") + 
  geom_vline(xintercept = mean(dist_0.9$icc), linetype = "dashed") + 
  geom_vline(xintercept = median(dist_0.9$icc), linetype = "dotted") + 
  labs(title = "Histogram of ICC Values for rho = .9 and cluster size 'U[30, 50]'" , x = "ICC Estimate", y = "Frequency") +
  theme_minimal()

7 Distributions of the Pooled ICC estimates

7.1 \(\rho = 0.05\):

results_all_converged |> filter(var_combo_graph == 0.05) |> 
  group_by(variance, model) |> 
  summarise(K = n(),
            "Mean Pooled Estimate" = mean(pooled_icc_est),
            var_est = var(pooled_icc_est),
            "Raw Bias Pooled Estimate" = mean(pooled_icc_est) - 0.05,
            "MCSE Raw Bias" = sqrt(var_est / K),
            RMSE = sqrt(sum((pooled_icc_est - 0.05)^2) / K),
            min_est = min(pooled_icc_est), 
            max_est = max(pooled_icc_est),
            .groups = "drop")|>
  kbl() |>
  kable_styling() |> 
  scroll_box(width = "900px", height = "500px")

variance	model	K	Mean Pooled Estimate	var_est	Raw Bias Pooled Estimate	MCSE Raw Bias	RMSE	min_est	max_est
Donner	REML	48000	0.0444794	1.93e-05	-0.0055206	2.00e-05	0.0070531	0.0194512	0.0689376
Donner	RVE	48000	0.0444223	1.95e-05	-0.0055777	2.02e-05	0.0071169	0.0197097	0.0696369
Fisher	REML	48000	0.0444866	1.93e-05	-0.0055134	2.00e-05	0.0070468	0.0197216	0.0690194
Fisher	RVE	48000	0.0444274	1.95e-05	-0.0055726	2.02e-05	0.0071105	0.0198642	0.0696210
Fisher TF	REML	48000	0.0461986	1.67e-05	-0.0038014	1.87e-05	0.0055821	0.0241912	0.0702339
Fisher TF	RVE	48000	0.0461986	1.67e-05	-0.0038014	1.87e-05	0.0055822	0.0241917	0.0702424
Hedges	REML	48000	0.0445867	1.92e-05	-0.0054133	2.00e-05	0.0069609	0.0196622	0.0691946
Hedges	RVE	48000	0.0445332	1.94e-05	-0.0054668	2.01e-05	0.0070212	0.0199373	0.0698602
Smith	REML	48000	0.0443622	1.94e-05	-0.0056378	2.01e-05	0.0071521	0.0195930	0.0687574
Smith	RVE	48000	0.0443009	1.96e-05	-0.0056991	2.02e-05	0.0072169	0.0196059	0.0693541
Swiger	REML	48000	0.0444893	1.93e-05	-0.0055107	2.00e-05	0.0070442	0.0197166	0.0690243
Swiger	RVE	48000	0.0444302	1.95e-05	-0.0055698	2.02e-05	0.0071078	0.0198623	0.0696227

results_all_converged |> filter(var_combo_graph == 0.05) |> 
#  group_by(var_combo_graph) |> 
  ggplot(aes(x = pooled_icc_est)) +
  geom_histogram( fill = "#404080", color = "#e9ecef", alpha = 0.6, position = 'identity', bins = 50) +
    facet_grid(model ~ variance ) 

7.2 \(\rho = 0.10\):

results_all_converged |> filter(var_combo_graph == 0.10) |> 
  group_by(variance, model) |> 
  summarise(K = n(),
            "Mean Pooled Estimate" = mean(pooled_icc_est),
            var_est = var(pooled_icc_est),
            "Raw Bias Pooled Estimate" = mean(pooled_icc_est) - 0.10,
            "MCSE Raw Bias" = sqrt(var_est / K),
            RMSE = sqrt(sum((pooled_icc_est - 0.10)^2) / K),
            min_est = min(pooled_icc_est), 
            max_est = max(pooled_icc_est), 
            .groups = "drop")|> 
  kbl() |>
  kable_styling() |> 
  scroll_box(width = "900px", height = "500px")

variance	model	K	Mean Pooled Estimate	var_est	Raw Bias Pooled Estimate	MCSE Raw Bias	RMSE	min_est	max_est
Donner	REML	48000	0.0929020	2.93e-05	-0.0070980	2.47e-05	0.0089287	0.0638537	0.1205014
Donner	RVE	48000	0.0929014	2.93e-05	-0.0070986	2.47e-05	0.0089284	0.0637351	0.1203858
Fisher	REML	48000	0.0928594	2.95e-05	-0.0071406	2.48e-05	0.0089687	0.0636333	0.1204919
Fisher	RVE	48000	0.0928520	2.95e-05	-0.0071480	2.48e-05	0.0089749	0.0635221	0.1203793
Fisher TF	REML	48000	0.0950363	2.59e-05	-0.0049637	2.32e-05	0.0071069	0.0673182	0.1239894
Fisher TF	RVE	48000	0.0950365	2.59e-05	-0.0049635	2.32e-05	0.0071069	0.0673183	0.1239837
Hedges	REML	48000	0.0930215	2.93e-05	-0.0069785	2.47e-05	0.0088291	0.0640675	0.1205711
Hedges	RVE	48000	0.0930298	2.92e-05	-0.0069702	2.47e-05	0.0088208	0.0639626	0.1204625
Smith	REML	48000	0.0926788	2.97e-05	-0.0073212	2.49e-05	0.0091247	0.0636260	0.1202392
Smith	RVE	48000	0.0926590	2.97e-05	-0.0073410	2.49e-05	0.0091417	0.0635153	0.1200836
Swiger	REML	48000	0.0928633	2.94e-05	-0.0071367	2.48e-05	0.0089651	0.0636440	0.1204934
Swiger	RVE	48000	0.0928562	2.94e-05	-0.0071438	2.48e-05	0.0089710	0.0635328	0.1203808

results_all_converged |> filter(var_combo_graph == 0.1) |> 
#  group_by(var_combo_graph) |> 
  ggplot(aes(x = pooled_icc_est)) +
  geom_histogram( fill = "#404080", color = "#e9ecef", alpha = 0.6, position = 'identity', bins = 50) +
    facet_grid(model ~ variance )

7.3 \(\rho = 0.15\):

results_all_converged |> filter(var_combo_graph == 0.15) |> 
  group_by(variance, model) |> 
  summarise(K = n(),
            "Mean Pooled Estimate" = mean(pooled_icc_est),
            var_est = var(pooled_icc_est),
            "Raw Bias Pooled Estimate" = mean(pooled_icc_est) - 0.15,
            "MCSE Raw Bias" = sqrt(var_est / K),
            RMSE = sqrt(sum((pooled_icc_est - 0.15)^2) / K),
            min_est = min(pooled_icc_est), 
            max_est = max(pooled_icc_est),
            .groups = "drop")|> 
  kbl() |>
  kable_styling() |> 
  scroll_box(width = "900px", height = "500px")

variance	model	K	Mean Pooled Estimate	var_est	Raw Bias Pooled Estimate	MCSE Raw Bias	RMSE	min_est	max_est
Donner	REML	48000	0.1418869	4.15e-05	-0.0081131	2.94e-05	0.0103601	0.1064030	0.1766836
Donner	RVE	48000	0.1418123	4.17e-05	-0.0081877	2.95e-05	0.0104256	0.1063915	0.1766836
Fisher	REML	48000	0.1418120	4.17e-05	-0.0081880	2.95e-05	0.0104287	0.1063188	0.1766414
Fisher	RVE	48000	0.1417274	4.19e-05	-0.0082726	2.95e-05	0.0105041	0.1063003	0.1766414
Fisher TF	REML	48000	0.1436102	3.84e-05	-0.0063898	2.83e-05	0.0089018	0.1106306	0.1773362
Fisher TF	RVE	48000	0.1436106	3.84e-05	-0.0063894	2.83e-05	0.0089014	0.1106319	0.1773445
Hedges	REML	48000	0.1420110	4.14e-05	-0.0079890	2.94e-05	0.0102592	0.1068142	0.1768486
Hedges	RVE	48000	0.1419499	4.15e-05	-0.0080501	2.94e-05	0.0103125	0.1068147	0.1768486
Smith	REML	48000	0.1416008	4.20e-05	-0.0083992	2.96e-05	0.0106104	0.1056585	0.1765680
Smith	RVE	48000	0.1414930	4.23e-05	-0.0085070	2.97e-05	0.0107085	0.1056643	0.1765680
Swiger	REML	48000	0.1418164	4.17e-05	-0.0081836	2.95e-05	0.0104245	0.1063331	0.1766487
Swiger	RVE	48000	0.1417324	4.19e-05	-0.0082676	2.95e-05	0.0104994	0.1063149	0.1766487

results_all_converged |> filter(var_combo_graph == 0.15) |> 
#  group_by(var_combo_graph) |> 
  ggplot(aes(x = pooled_icc_est)) +
  geom_histogram( fill = "#404080", color = "#e9ecef", alpha = 0.6, position = 'identity', bins = 50) +
    facet_grid(model~ variance )

7.4 \(\rho = 0.25\):

results_all_converged |> filter(var_combo_graph == 0.25) |> 
  group_by(variance, model) |> 
  summarise(K = n(),
            "Mean Pooled Estimate" = mean(pooled_icc_est),
            var_est = var(pooled_icc_est),
            "Raw Bias Pooled Estimate" = mean(pooled_icc_est) - 0.25,
            "MCSE Raw Bias" = sqrt(var_est / K),
            RMSE = sqrt(sum((pooled_icc_est - 0.25)^2) / K),
            min_est = min(pooled_icc_est), 
            max_est = max(pooled_icc_est),
            .groups = "drop")|>

  kbl() |>
  kable_styling() |> 
  scroll_box(width = "900px", height = "500px")

variance	model	K	Mean Pooled Estimate	var_est	Raw Bias Pooled Estimate	MCSE Raw Bias	RMSE	min_est	max_est
Donner	REML	48000	0.2413431	6.46e-05	-0.0086569	3.67e-05	0.0118127	0.2013939	0.2871551
Donner	RVE	48000	0.2411995	6.50e-05	-0.0088005	3.68e-05	0.0119352	0.2016610	0.2870117
Fisher	REML	48000	0.2412389	6.49e-05	-0.0087611	3.68e-05	0.0119028	0.2012907	0.2869823
Fisher	RVE	48000	0.2410835	6.54e-05	-0.0089165	3.69e-05	0.0120371	0.2013082	0.2868213
Fisher TF	REML	48000	0.2410603	6.57e-05	-0.0089397	3.70e-05	0.0120662	0.1997265	0.2822733
Fisher TF	RVE	48000	0.2410612	6.57e-05	-0.0089388	3.70e-05	0.0120653	0.1997266	0.2822763
Hedges	REML	48000	0.2414566	6.45e-05	-0.0085434	3.67e-05	0.0117264	0.2011547	0.2875039
Hedges	RVE	48000	0.2413259	6.49e-05	-0.0086741	3.68e-05	0.0118374	0.2017192	0.2873750
Smith	REML	48000	0.2410222	6.54e-05	-0.0089778	3.69e-05	0.0120820	0.2012777	0.2869650
Smith	RVE	48000	0.2408420	6.60e-05	-0.0091580	3.71e-05	0.0122404	0.2012916	0.2868021
Swiger	REML	48000	0.2412436	6.49e-05	-0.0087564	3.68e-05	0.0118986	0.2012878	0.2869940
Swiger	RVE	48000	0.2410888	6.54e-05	-0.0089112	3.69e-05	0.0120324	0.2013215	0.2868341

results_all_converged |> filter(var_combo_graph == 0.25) |> 
#  group_by(var_combo_graph) |> 
  ggplot(aes(x = pooled_icc_est)) +
  geom_histogram( fill = "#404080", color = "#e9ecef", alpha = 0.6, position = 'identity', bins = 50) +
    facet_grid(model~ variance )

7.5 \(\rho = 0.50\):

results_all_converged |> filter(var_combo_graph == 0.50) |> 
  group_by(variance, model) |> 
  summarise(K = n(),
            "Mean Pooled Estimate" = mean(pooled_icc_est),
      #      "Median Pooled Estimate" = median(pooled_icc_est),
            var_est = var(pooled_icc_est),
            "Raw Bias Pooled Estimate" = mean(pooled_icc_est) - 0.50,
            "MCSE Raw Bias" = sqrt(var_est / K),
            RMSE = sqrt(sum((pooled_icc_est - 0.50)^2) / K), 
            min_est = min(pooled_icc_est), 
            max_est = max(pooled_icc_est),
            .groups = "drop")|>
  kbl() |>
  kable_styling() |> 
  scroll_box(width = "900px", height = "500px")

variance	model	K	Mean Pooled Estimate	var_est	Raw Bias Pooled Estimate	MCSE Raw Bias	RMSE	min_est	max_est
Donner	REML	48000	0.4953113	0.0000912	-0.0046887	4.36e-05	0.0106376	0.4411449	0.5419561
Donner	RVE	48000	0.4953182	0.0000913	-0.0046818	4.36e-05	0.0106394	0.4411045	0.5420179
Fisher	REML	48000	0.4952857	0.0000913	-0.0047143	4.36e-05	0.0106546	0.4410808	0.5419028
Fisher	RVE	48000	0.4952923	0.0000914	-0.0047077	4.36e-05	0.0106570	0.4410339	0.5419671
Fisher TF	REML	48000	0.4880025	0.0001072	-0.0119975	4.73e-05	0.0158473	0.4341028	0.5367767
Fisher TF	RVE	48000	0.4880040	0.0001072	-0.0119960	4.73e-05	0.0158458	0.4340954	0.5367603
Hedges	REML	48000	0.4953049	0.0000911	-0.0046951	4.36e-05	0.0106353	0.4412048	0.5418731
Hedges	RVE	48000	0.4953117	0.0000912	-0.0046883	4.36e-05	0.0106370	0.4411663	0.5419321
Smith	REML	48000	0.4952682	0.0000915	-0.0047318	4.37e-05	0.0106725	0.4411312	0.5421544
Smith	RVE	48000	0.4952744	0.0000916	-0.0047256	4.37e-05	0.0106756	0.4409836	0.5422489
Swiger	REML	48000	0.4952861	0.0000913	-0.0047139	4.36e-05	0.0106541	0.4410833	0.5419024
Swiger	RVE	48000	0.4952927	0.0000914	-0.0047073	4.36e-05	0.0106565	0.4410366	0.5419665

results_all_converged |> filter(var_combo_graph == 0.5) |> 
#  group_by(var_combo_graph) |> 
  ggplot(aes(x = pooled_icc_est)) +
  geom_histogram( fill = "#404080", color = "#e9ecef", alpha = 0.6, position = 'identity', bins = 50) +
    facet_grid(model ~ variance )

7.6 \(\rho = 0.90\):

results_all_converged |> filter(var_combo_graph == 0.90) |> 
  group_by(variance, model) |> 
  summarise(K = n(),
            "Mean Pooled Estimate" = mean(pooled_icc_est),
      #      "Median Pooled Estimate" = median(pooled_icc_est),
            var_est = var(pooled_icc_est),
            "Raw Bias Pooled Estimate" = mean(pooled_icc_est) - 0.90,
            "MCSE Raw Bias" = sqrt(var_est / K),
            RMSE = sqrt(sum((pooled_icc_est - 0.90)^2) / K), 
            min_est = min(pooled_icc_est), 
            max_est = max(pooled_icc_est),
            .groups = "drop")|>
  kbl() |>
  kable_styling() |> 
  scroll_box(width = "900px", height = "500px")

variance	model	K	Mean Pooled Estimate	var_est	Raw Bias Pooled Estimate	MCSE Raw Bias	RMSE	min_est	max_est
Donner	REML	48000	0.9026873	1.98e-05	0.0026873	2.03e-05	0.0051952	0.8735445	0.9278472
Donner	RVE	48000	0.9026360	1.98e-05	0.0026360	2.03e-05	0.0051716	0.8732500	0.9277722
Fisher	REML	48000	0.9027883	1.99e-05	0.0027883	2.03e-05	0.0052565	0.8735703	0.9280164
Fisher	RVE	48000	0.9027390	1.99e-05	0.0027390	2.03e-05	0.0052311	0.8732787	0.9279004
Fisher TF	REML	48000	0.8970148	2.40e-05	-0.0029852	2.24e-05	0.0057375	0.8673503	0.9201010
Fisher TF	RVE	48000	0.8970153	2.40e-05	-0.0029847	2.24e-05	0.0057371	0.8673525	0.9201015
Hedges	REML	48000	0.9025850	1.98e-05	0.0025850	2.03e-05	0.0051422	0.8735030	0.9274451
Hedges	RVE	48000	0.9025347	1.98e-05	0.0025347	2.03e-05	0.0051211	0.8732056	0.9274583
Smith	REML	48000	0.9029853	2.00e-05	0.0029853	2.04e-05	0.0053753	0.8738152	0.9280318
Smith	RVE	48000	0.9029397	2.00e-05	0.0029397	2.04e-05	0.0053475	0.8735355	0.9279116
Swiger	REML	48000	0.9027842	1.99e-05	0.0027842	2.03e-05	0.0052539	0.8735697	0.9280039
Swiger	RVE	48000	0.9027347	1.99e-05	0.0027347	2.03e-05	0.0052285	0.8732780	0.9278907

results_all_converged |> filter(var_combo_graph == 0.9) |> 
#  group_by(var_combo_graph) |> 
  ggplot(aes(x = pooled_icc_est)) +
  geom_histogram( fill = "#404080", color = "#e9ecef", alpha = 0.6, position = 'identity', bins = 50) +
    facet_grid(model ~ variance )

7.7 Skewness and Kurtosis

perf_crit |> 
  group_by(variance, model, var_combo_graph) |> 
  summarise(m_s = mean(skewness),
          #  min_s = min(skewness),
          #  max_s = max(skewness),
            m_k = mean(kurtosis),
          #  min_k = min(kurtosis),
          #  max_k = max(kurtosis)
          .groups = "drop")|>
  pivot_wider(names_from = variance, values_from = c(m_s, m_k)) |>
  kable(booktabs = TRUE, col.names = c("Model",
                    "ICC Value",
                    "Donner",
                    "Fisher",
                    "Fisher TF",
                    "Hedges",
                    "Smith",
                    "Swiger",
                    "Donner",
                    "Fisher",
                    "Fisher TF",
                    "Hedges",
                    "Smith",
                    "Swiger"
                    )) |>
  kable_styling() |> 
  add_header_above(c(" " = 2, "Mean Skewness" = 6, "Mean Kurtosis" = 6)) |> 
  column_spec(9:14, bold = T, background = "grey") 

		Mean Skewness						Mean Kurtosis
Model	ICC Value	Donner	Fisher	Fisher TF	Hedges	Smith	Swiger	Donner	Fisher	Fisher TF	Hedges	Smith	Swiger
REML	0.05	0.0592022	0.0594254	0.0866294	0.0593649	0.0591826	0.0594198	2.954853	2.955471	2.953104	2.954562	2.955909	2.955497
REML	0.10	0.0723233	0.0724675	0.0749330	0.0718898	0.0727667	0.0724434	3.001821	3.001935	2.995341	3.001267	3.002176	3.001900
REML	0.15	0.0357151	0.0359730	0.0377721	0.0360568	0.0364047	0.0359696	2.947726	2.948172	2.955624	2.947198	2.948294	2.948137
REML	0.25	0.0323212	0.0323009	0.0280801	0.0319869	0.0327775	0.0322975	3.010076	3.009571	3.014416	3.008855	3.010905	3.009555
REML	0.50	-0.0362609	-0.0361881	-0.0318783	-0.0357684	-0.0362934	-0.0361812	2.980498	2.980601	2.978749	2.980647	2.980997	2.980596
REML	0.90	-0.0973374	-0.0970725	-0.1078941	-0.0969041	-0.0974283	-0.0970904	3.027450	3.027156	3.022494	3.027109	3.026623	3.027172
RVE	0.05	0.0684156	0.0684702	0.0866546	0.0682543	0.0689004	0.0684497	2.951688	2.952241	2.953011	2.951607	2.952394	2.952227
RVE	0.10	0.0742096	0.0746966	0.0749193	0.0736681	0.0752857	0.0746794	2.998217	2.998166	2.995317	2.997296	2.998940	2.998164
RVE	0.15	0.0359919	0.0364623	0.0378290	0.0359824	0.0374692	0.0364526	2.949006	2.949565	2.955679	2.948395	2.950002	2.949530
RVE	0.25	0.0336099	0.0338578	0.0281223	0.0329711	0.0346274	0.0338303	3.005622	3.004860	3.014341	3.005028	3.005913	3.004849
RVE	0.50	-0.0364144	-0.0363523	-0.0318745	-0.0358920	-0.0365080	-0.0363443	2.980415	2.980558	2.978801	2.980482	2.981113	2.980552
RVE	0.90	-0.1068451	-0.1071497	-0.1078788	-0.1061320	-0.1078974	-0.1071320	3.024866	3.024611	3.022530	3.024555	3.023916	3.024628

# %>%
#   row_spec(3:5, bold = T, color = "white", background = "#D7261E")


perf_crit |> 
  filter(var_combo_graph == 0.05) |>
  group_by(model, variance) |> 
  summarise(
            m_b = mean(bias_icc),
            m_rb = mean(rel_bias_icc),
            m_s = mean(skewness),
          #  min_s = min(skewness),
          #  max_s = max(skewness),
            m_k = mean(kurtosis),
          #  min_k = min(kurtosis),
          #  max_k = max(kurtosis)
          .groups = "drop")|>
  kable(booktabs = TRUE, col.names = c("Model",
                    "ICC Value", "Mean Bias", "Mean Relative Bias", "Mean Skewness", "Mean Kurtosis")) |>
  kable_styling() |> 
  add_header_above(c(" " = 2, "ICC = 0.05" = 4))

		ICC = 0.05
Model	ICC Value	Mean Bias	Mean Relative Bias	Mean Skewness	Mean Kurtosis
REML	Donner	-0.0055206	-0.1104120	0.0592022	2.954853
REML	Fisher	-0.0055134	-0.1102687	0.0594254	2.955471
REML	Fisher TF	-0.003801 TF	-0.003801

perf_crit |> 
  filter(var_combo_graph == 0.1) |>
  group_by(model, variance) |> 
  summarise(
            m_b = mean(bias_icc),
            m_rb = mean(rel_bias_icc),
            m_s = mean(skewness),
          #  min_s = min(skewness),
          #  max_s = max(skewness),
            m_k = mean(kurtosis),
          #  min_k = min(kurtosis),
          #  max_k = max(kurtosis)
          .groups = "drop")|>
  kable(booktabs = TRUE, col.names = c("Model",
                    "ICC Value", "Mean Bias", "Mean Relative Bias", "Mean Skewness", "Mean Kurtosis")) |>
  kable_styling() |> 
  add_header_above(c(" " = 2, "ICC = 0.10" = 4))

		ICC = 0.10
Model	ICC Value	Mean Bias	Mean Relative Bias	Mean Skewness	Mean Kurtosis
REML	Donner	-0.0070980	-0.0709803	0.0723233	3.001821
REML	Fisher	-0.0071406	-0.0714056	0.0724675	3.001935
REML	Fisher TF	-0.004963 TF	-0.004963

perf_crit |> 
  filter(var_combo_graph == 0.15) |>
  group_by(model, variance) |> 
  summarise(
            m_b = mean(bias_icc),
            m_rb = mean(rel_bias_icc),
            m_s = mean(skewness),
          #  min_s = min(skewness),
          #  max_s = max(skewness),
            m_k = mean(kurtosis),
          #  min_k = min(kurtosis),
          #  max_k = max(kurtosis)
          .groups = "drop")|>
  kable(booktabs = TRUE, col.names = c("Model",
                    "ICC Value", "Mean Bias", "Mean Relative Bias", "Mean Skewness", "Mean Kurtosis")) |>
  kable_styling() |> 
  add_header_above(c(" " = 2, "ICC = 0.15" = 4))

		ICC = 0.15
Model	ICC Value	Mean Bias	Mean Relative Bias	Mean Skewness	Mean Kurtosis
REML	Donner	-0.0081131	-0.0540873	0.0357151	2.947726
REML	Fisher	-0.0081880	-0.0545870	0.0359730	2.948172
REML	Fisher TF	-0.006389 TF	-0.006389

perf_crit |> 
  filter(var_combo_graph == 0.25) |>
  group_by(model, variance) |> 
  summarise(
            m_b = mean(bias_icc),
            m_rb = mean(rel_bias_icc),
            m_s = mean(skewness),
          #  min_s = min(skewness),
          #  max_s = max(skewness),
            m_k = mean(kurtosis),
          #  min_k = min(kurtosis),
          #  max_k = max(kurtosis)
          .groups = "drop")|>
  kable(booktabs = TRUE, col.names = c("Model",
                    "ICC Value", "Mean Bias", "Mean Relative Bias", "Mean Skewness", "Mean Kurtosis")) |>
  kable_styling() |> 
  add_header_above(c(" " = 2, "ICC = 0.25" = 4))

		ICC = 0.25
Model	ICC Value	Mean Bias	Mean Relative Bias	Mean Skewness	Mean Kurtosis
REML	Donner	-0.0086569	-0.0346275	0.0323212	3.010076
REML	Fisher	-0.0087611	-0.0350445	0.0323009	3.009571
REML	Fisher TF	-0.008939 TF	-0.008938

perf_crit |> 
  filter(var_combo_graph == 0.5) |>
  group_by(model, variance) |> 
  summarise(
            m_b = mean(bias_icc),
            m_rb = mean(rel_bias_icc),
            m_s = mean(skewness),
          #  min_s = min(skewness),
          #  max_s = max(skewness),
            m_k = mean(kurtosis),
          #  min_k = min(kurtosis),
          #  max_k = max(kurtosis)
          .groups = "drop")|>
  kable(booktabs = TRUE, col.names = c("Model",
                    "ICC Value", "Mean Bias", "Mean Relative Bias", "Mean Skewness", "Mean Kurtosis")) |>
  kable_styling() |> 
  add_header_above(c(" " = 2, "ICC = 0.50" = 4))

		ICC = 0.50
Model	ICC Value	Mean Bias	Mean Relative Bias	Mean Skewness	Mean Kurtosis
REML	Donner	-0.0046887	-0.0093774	-0.0362609	2.980498
REML	Fisher	-0.0047143	-0.0094286	-0.0361881	2.980601
REML	Fisher TF	-0.011997 TF	-0.011996

perf_crit |> 
  filter(var_combo_graph == 0.9) |>
  group_by(model, variance) |> 
  summarise(
            m_b = mean(bias_icc),
            m_rb = mean(rel_bias_icc),
            m_s = mean(skewness),
          #  min_s = min(skewness),
          #  max_s = max(skewness),
            m_k = mean(kurtosis),
          #  min_k = min(kurtosis),
          #  max_k = max(kurtosis)
          .groups = "drop")|>
  kable(booktabs = TRUE, col.names = c("Model",
                    "ICC Value", "Mean Bias", "Mean Relative Bias", "Mean Skewness", "Mean Kurtosis")) |>
  kable_styling() |> 
  add_header_above(c(" " = 2, "ICC = 0.90" = 4))

		ICC = 0.90
Model	ICC Value	Mean Bias	Mean Relative Bias	Mean Skewness	Mean Kurtosis
REML	Donner	0.0026873	0.0029859	-0.0973374	3.027450
REML	Fisher	0.0027883	0.0030982	-0.0970725	3.027156
REML	Fisher TF	-0.002985 TF	-0.002984

perf_crit_agg <- perf_crit |> 
  group_by(model, variance, icc_est_n , var_combo_graph, tau) |> 
  summarise(skewness_mean = mean(skewness), 
            kurtosis_mean = mean(kurtosis),
            n = n(),
            .groups = "drop")



ID1 <- c(0.05, 0.10, 0.15)

perf_crit_agg |> filter(var_combo_graph < 0.2) |> ggplot(aes( var_combo_graph, skewness_mean, col=variance, shape = model ) ) + 
  facet_grid(tau ~ icc_est_n, 
             labeller = label_bquote(rows = tau == .(tau), 
                                            cols = k ==.(icc_est_n))) +

  geom_point( alpha=0.75 ) + 
  geom_line( alpha=0.75 ) +
  geom_hline( yintercept = 0 ) + 
  scale_x_continuous("ICC Parameter", labels = as.character(ID1), breaks = ID1) +
  theme_minimal()

ID2 <- c(0.25, 0.5, 0.9)

perf_crit_agg |> filter(var_combo_graph > 0.2) |> ggplot(aes( var_combo_graph, skewness_mean, col=variance, shape = model ) ) + 
  facet_grid(tau ~ icc_est_n, 
             labeller = label_bquote(rows = tau == .(tau), 
                                            cols = k ==.(icc_est_n))) +

  geom_point( alpha=0.75 ) + 
  geom_line( alpha=0.75 ) +
  geom_hline( yintercept = 0 ) + 
  scale_x_continuous("ICC Parameter", labels = as.character(ID2), breaks = ID2) +
  theme_minimal()

ID1 <- c(0.05, 0.10, 0.15)

perf_crit_agg |> filter(var_combo_graph < 0.2) |> ggplot(aes( var_combo_graph, kurtosis_mean, col=variance, shape = model ) ) + 
  facet_grid(tau ~ icc_est_n, 
             labeller = label_bquote(rows = tau == .(tau), 
                                            cols = k ==.(icc_est_n))) +

  geom_point( alpha=0.75 ) + 
  geom_line( alpha=0.75 ) +
  geom_hline( yintercept = 3 ) + 
  scale_x_continuous("ICC Parameter", labels = as.character(ID1), breaks = ID1) +
  theme_minimal()

ID2 <- c(0.25, 0.5, 0.9)

perf_crit_agg |> filter(var_combo_graph > 0.2) |> ggplot(aes( var_combo_graph, kurtosis_mean, col=variance, shape = model ) ) + 
  facet_grid(tau ~ icc_est_n, 
             labeller = label_bquote(rows = tau == .(tau), 
                                            cols = k ==.(icc_est_n))) +

  geom_point( alpha=0.75 ) + 
  geom_line( alpha=0.75 ) +
  geom_hline( yintercept = 3 ) + 
  scale_x_continuous("ICC Parameter", labels = as.character(ID2), breaks = ID2) +
  theme_minimal()

Other Graphs

perf_crit |> 
  ggplot(aes(x = factor(var_combo_graph), y = skewness, color = factor(icc_est_n))) +
  geom_point(alpha = .5, position = position_jitter(width = .3)) + 
  geom_hline(yintercept = 0, linetype = "solid")  + 
  facet_grid(variance ~ model, 
             labeller = label_bquote(rows = .(variance), 
                                            cols = .(model))) +
  labs(x = "ICC Value", y = "Skewness", color = "k") 

perf_crit |> 
  ggplot(aes(x = factor(var_combo_graph), y = skewness, color = factor(tau))) +
  geom_point(alpha = .5, position = position_jitter(width = .3)) +
    geom_hline(yintercept = 0, linetype = "solid")  + 
  facet_grid(variance ~ icc_est_n, 
             labeller = label_bquote(rows = .(variance), 
                                            cols = k ==.(icc_est_n))) +
  labs(x = "ICC Value", y = "Skewness", color = "tau") 

perf_crit |> 
  ggplot(aes(x = factor(var_combo_graph), y = skewness, color = nj_graph)) +
  geom_point(alpha = .5, position = position_jitter(width = .3)) +
    geom_hline(yintercept = 0, linetype = "solid")  + 
  facet_grid(variance ~ icc_est_n, 
             labeller = label_bquote(rows = .(variance), 
                                            cols = k ==.(icc_est_n))) +
  labs(x = "ICC Value", y = "Skewness", color = "") 

perf_crit |> 
  ggplot(aes(x = factor(var_combo_graph), y = skewness, color = n_bar_graph)) +
  geom_point(alpha = .5, position = position_jitter(width = .3)) +
  geom_hline(yintercept = 0, linetype = "solid")  + 
  facet_grid(variance ~ icc_est_n, 
             labeller = label_bquote(rows = .(variance), 
                                            cols = k ==.(icc_est_n))) +
  labs(x = "ICC Value", y = "Skewness", color = "") 

perf_crit |> 
  ggplot(aes(x = factor(var_combo_graph), y = kurtosis, color = factor(icc_est_n))) +
  geom_point(alpha = .5, position = position_jitter(width = .3)) +
  geom_hline(yintercept = 3, linetype = "solid")  + 
  facet_grid(variance ~ model) + 
  labs(x = "ICC Value", y = "Kurtosis", color = "k") 

perf_crit |> 
  ggplot(aes(x = factor(var_combo_graph), y = kurtosis, color = factor(tau))) +
  geom_point(alpha = .5, position = position_jitter(width = .3)) +
    geom_hline(yintercept = 3, linetype = "solid")  + 
  facet_grid(variance ~ icc_est_n, 
             labeller = label_bquote(rows = .(variance), 
                                            cols = k ==.(icc_est_n))) +
  labs(x = "ICC Value", y = "Kurtosis", color = "tau") 

perf_crit |> 
  ggplot(aes(x = factor(var_combo_graph), y = kurtosis, color = nj_graph)) +
  geom_point(alpha = .5, position = position_jitter(width = .3)) +
    geom_hline(yintercept = 3, linetype = "solid")  + 
  facet_grid(variance ~ icc_est_n, 
             labeller = label_bquote(rows = .(variance), 
                                            cols = k ==.(icc_est_n))) +
  labs(x = "ICC Value", y = "Kurtosis", color = "") 

perf_crit |> 
  ggplot(aes(x = factor(var_combo_graph), y = kurtosis, color = n_bar_graph)) +
  geom_point(alpha = .5, position = position_jitter(width = .3)) +
  geom_hline(yintercept = 3, linetype = "solid")  + 
  facet_grid(variance ~ icc_est_n, 
             labeller = label_bquote(rows = .(variance), 
                                            cols = k ==.(icc_est_n))) +
  labs(x = "ICC Value", y = "Kurtosis", color = "")